//Environmental Regulation and Product Attributes: The Case of European Passenger Vehicle Greenhouse Gas Emissions Standards
//Yujie Lin, Joshua Linn
//Stata/MP 17.0
//Apr 25, 2022

* This file generates all estimation results and figures for the main text and appendix in the paper.


* Notes: Put the data and the following codes in the same folder before running the codes.


**********************************************************************************************************************************
**#                        Part I: Main Text
**********************************************************************************************************************************
********************************************************************************
**# Table 1: Summary Statistic
********************************************************************************
use "data.dta", clear
sum price
reghdfe logshare_jcy price tax energycost1 loghw logw logsize  ,   ///
	absorb(i.fuelcat_harmonized2 i.groupmodel2#i.bodytype_harmonized4 i.country2#i.year i.segment2#i.year)  vce(cluster group_mdltrim)
est store ols
gen sample=e(sample)
tab post if sample==1

tabstat reg co2emit price tax enginehorsepower grossvehicleweight size energycost1 fuelcons enginecyl if sample==1,by(post) stat(mean sd)

//market share of plug-in vehicles in each time period
tab fuelcat_harmonized2 post 

drop reg_cy  
bysort country year: egen reg_cy=total(reg), missing
bysort country year fuelcat_harmonized2: egen reg_cyf=total(reg), missing
bysort country year fuelcat_harmonized2: gen share_cyf=reg_cyf/reg_cy

bysort year: egen reg_y=total(reg), missing
bysort year fuelcat_harmonized2: egen reg_yf=total(reg), missing
bysort year fuelcat_harmonized2: gen share_yf=reg_yf/reg_y

bysort post: egen reg_p=total(reg), missing
bysort post fuelcat_harmonized2: egen reg_pf=total(reg), missing
bysort post fuelcat_harmonized2: gen share_pf=reg_pf/reg_p

sum share_cyf share_yf share_pf if fuelcat_harmonized2==2 & post==1 
sum share_cyf share_yf share_pf if fuelcat_harmonized2==2 & post==2
sum share_cyf share_yf share_pf if fuelcat_harmonized2==2 & post==3 



********************************************************************************
**# Figure 1: co2 trend
********************************************************************************
tabstat co2emit [aw=reg], by(year) stats(p25 p50 p75)
//save the data to a new data file and we can draw the figure
/*
    year |       p25       p50       p75
---------+------------------------------
    2005 |       135     150.5       166
    2006 |       133       149  164.9999
    2007 |   130.455  146.3699  162.1424
    2008 |  125.8095  142.1889  158.4926
    2009 |       120       138  152.8889
    2010 |  118.1646       130  145.8485
    2011 |  110.3419  124.2873    138.22
    2012 |       109  119.6923       134
    2013 |       104  116.8261       128
    2014 |       102       114       124
    2015 |        99       109       120
    2016 |        99       108       119
    2017 |        99       108       119
---------+------------------------------
   Total |       113       128     148.4
----------------------------------------
*/

********************************************************************************
**# Table 2: demand estimation
********************************************************************************
use "data.dta", clear

//the 1% upper side outliers of price
sum price [fw=reg], d
gen routlier3_price=cond(price>r(p99),1,0)
sum price [fw=reg], d
gen routlier4_price=cond(price>r(p95),1,0)
drop if routlier3_price==1 

sum enginehorsepower grossvehicleweight loghw loghw2
drop loghw
gen loghw=log(enginehorsepower*100/(grossvehicleweight*1000))
label var loghw "log(enginehorsepower*100/(grossvehicleweight*1000))"

* 1.1 MNL_OLS
drop fe_* gamma_* wtp_* resid_* dep_*
reghdfe logshare_jcy price tax energycost1 loghw logw logsize  ,   ///
	absorb(fe_cy_ols=i.country2#i.year fe_mb_ols=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_ols=i.segment2#i.year, resid(resid_ols))  vce(cluster group_mdltrim)
est store ols
matrix beta_ols=e(b)'
matrix list beta_ols    
predict fe_ols, d    

gen gamma_jcy_ols=logshare_jcy-beta_ols[1,1]*price-beta_ols[2,1]*tax-beta_ols[3,1]*energycost1    ///
	-beta_ols[4,1]*loghw-beta_ols[5,1]*logw-fe_cy_ols-beta_ols[7,1]
gen wtp_jcy_ols=gamma_jcy_ols/(-1*beta_ols[1,1])

sum gamma* wtp*

//individual elasticity
gen elastic_price_ols=beta_ols[1,1]*price
sum elastic_price_ols,d
//aggregated elasticity
bysort country year: asgen avgela_price_ols=elastic_price_ols, weight(reg)
sum avgela_price_ols,d

//WTP
matrix wtp_ols=beta_ols/abs(beta_ols[1,1])
matrix list wtp_ols


* 1.2 MNL_IV
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize    ///
	(price=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff),   ///
	absorb(fe_cy_mnliv=i.country2#i.year fe_mb_mnliv=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_mnliv=i.segment2#i.year, resid(resid_mnliv)) first ffirst cluster(group_mdltrim)
est store mnliv
matrix beta_mnliv=e(b)'
matrix list beta_mnliv

predict fe_mnliv, d

* add constant
matrix list beta_mnliv
gen dep_mnl=logshare_jcy-beta_mnliv[1,1]*price
reghdfe dep_mnl tax energycost1 loghw logw logsize  ,   ///
	absorb(fe_cy_mnliv2=i.country2#i.year fe_mb_mnliv2=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_mnliv2=i.segment2#i.year, resid(resid_mnliv2))  vce(cluster group_mdltrim)
est store mnliv2
matrix beta_mnliv2=e(b)'
matrix list beta_mnliv2
predict fe_mnliv2,d
sum fe_cy_mnliv  

* gamma and wtp
gen gamma_jcy_mnliv=logshare_jcy-beta_mnliv[1,1]*price-beta_mnliv[2,1]*tax-beta_mnliv[3,1]*energycost1    ///
	-beta_mnliv[4,1]*loghw-beta_mnliv[5,1]*logw-fe_cy_mnliv
gen wtp_jcy_mnliv=gamma_jcy_mnliv/beta_mnliv[1,1]*-1

gen gamma_jcy_mnliv2=dep_mnl-beta_mnliv2[1,1]*tax-beta_mnliv2[2,1]*energycost1-beta_mnliv2[3,1]*loghw    ///
	-beta_mnliv2[4,1]*logw-fe_cy_mnliv2-beta_mnliv2[6,1]
gen wtp_jcy_mnliv2=gamma_jcy_mnliv2/beta_mnliv[1,1]*-1

//elasticity and WTP
gen elastic_price_mnl=(1-share_jcy)*beta_mnliv[1,1]*price
sum elastic_price_mnl,d
bysort country year: asgen avgela_price_mnl=elastic_price_mnl, weight(reg)
sum avgela_price_mnl,d

matrix wtp_mnl=beta_mnliv/abs(beta_mnliv[1,1])
matrix list wtp_mnl


* 1.3 NL_IV
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_nl=i.country2#i.year fe_mb_nl=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_nl=i.segment2#i.year, resid(resid_nl)) first ffirst cluster(group_mdltrim)
est store nl
matrix beta=e(b)'
matrix list beta

predict fe_nl, d

* add constant
matrix list beta
gen dep_nl=logshare_jcy-beta[1,1]*price-beta[2,1]*logshare_cyjinso-beta[3,1]*logshare_cyoins
reghdfe dep_nl tax energycost1 loghw logw logsize  ,   ///
	absorb(fe_cy_nlt=i.country2#i.year fe_mb_nlt=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_nlt=i.segment2#i.year, resid(resid_nlt))  vce(cluster group_mdltrim)
est store nlt
matrix betanlt=e(b)'
matrix list betanlt
predict fe_nlt, d
sum fe_cy_nl   

* gamma and wtp
gen gamma_jcy1=logshare_jcy-beta[1,1]*price-beta[2,1]*logshare_cyjinso-beta[3,1]*logshare_cyoins    ///
	-beta[4,1]*tax-beta[5,1]*energycost1-beta[6,1]*loghw-beta[7,1]*logw-fe_cy_nl
gen wtp_jcy1=gamma_jcy1/beta[1,1]*-1

gen gamma_jcy1t=dep_nl-betanlt[1,1]*tax-betanlt[2,1]*energycost1-betanlt[3,1]*loghw-    ///
	betanlt[4,1]*logw-fe_cy_nlt-betanlt[6,1]
gen wtp_jcy1t=gamma_jcy1t/beta[1,1]*-1

sum gamma* wtp* 
sum logshare_jcy
sum fe_mb* fe_cy* resid_*

//elasticity and wtp for attributes
gen elastic_price=((1-share_cyjinso)/(1-beta[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta[1,1]*price
sum elastic_price,d
bysort country year: asgen avgela_price=elastic_price, weight(reg)
sum avgela_price,d	

matrix wtp=beta/abs(beta[1,1])
matrix list wtp

//quality index by country-year
bysort country year: egen gamma_cy_ols=mean(gamma_jcy_ols)
bysort country year: egen wtp_cy_ols=mean(wtp_jcy_ols)
bysort country year: egen gamma_cy_mnliv=mean(gamma_jcy_mnliv)
bysort country year: egen wtp_cy_mnliv=mean(wtp_jcy_mnliv)
bysort country year: egen gamma_cy_mnliv2=mean(gamma_jcy_mnliv2)
bysort country year: egen wtp_cy_mnliv2=mean(wtp_jcy_mnliv2)
bysort country year: egen gamma_cy1=mean(gamma_jcy1)
bysort country year: egen wtp_cy1=mean(wtp_jcy1)
bysort country year: egen gamma_cy1t=mean(gamma_jcy1t)
bysort country year: egen wtp_cy1t=mean(wtp_jcy1t)
bysort year: egen gamma_y1=mean(gamma_jcy1)
bysort year: egen wtp_y1=mean(wtp_jcy1)

*************************************
* output for Table 2
*************************************
esttab ols mnliv nl, b(%7.3f) se(%7.3f) stats(N r2_a) br compress nostar    ///
	order(price logshare_cyjinso logshare_cyoins tax energycost1 loghw logw logsize _cons)
estimates tab ols mnliv nl, b(%7.3f) se(%7.3f) stats(N r2_a) 

unique vehicle if _est_nl==1
unique groupmodel2 if _est_nl==1

*************************************
* output for Table 3: average price elasticity 
*************************************
sum avgela_price_ols, d 
matrix A=(r(p50)\r(mean)\r(sd)\r(p5)\r(p95))
foreach x in avgela_price_mnl avgela_price{
    sum `x', d
	matrix A=(A,(r(p50)\r(mean)\r(sd)\r(p5)\r(p95)))
}
matrix list A
matrix drop A



********************************************************************************
**# Table 4: DID regressions
********************************************************************************
use "data.dta", clear
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(i.country2#i.year i.groupmodel2#i.bodytype_harmonized4    ///
	 i.segment2#i.year) first ffirst cluster(group_mdltrim)
est store nl
matrix beta=e(b)'
matrix list beta
gen snl=e(sample)

//normalize the loghw and logw
drop loghw3 logw2
gen loghw3=loghw/abs(beta[1,1])
gen logw2=logw/abs(beta[1,1])
label var loghw3 "loghw (hp/kg) normalized by the absolute value of price coefficient"
label var logw2 "logw (ton) normalized by the absolute value of price coefficient"

* main regressions using firm level stringnecy
foreach x in wtp_jcy1 loghw3 logw2 logprice {
	reghdfe `x' c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3   ///
		[fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year) vce(cluster group_mdltrim) 
	est store wf0_`x'
	count if e(sample)==1
	matrix n_wf0_`x'=r(N)
	test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
	matrix test_wf0_`x'=(r(drop),  r(df_r), r(F), r(df), r(p))'
}

* output
esttab wf0_wtp_jcy1 wf0_loghw3 wf0_logw2 wf0_logprice,     ///
	b(%7.3f) se(%7.3f) stats(N r2_a) br compress order(*.post#* stringency_* _cons) nostar nogap

matrix A=(test_wf0_wtp_jcy1, test_wf0_loghw3, test_wf0_logw2, test_wf0_logprice)
matrix list A

matrix B=(n_wf0_wtp_jcy1, n_wf0_loghw3, n_wf0_logw2, n_wf0_logprice)
matrix list B

matrix drop A B






**********************************************************************************************************************************
**#                        Part II: Appendix-----Demand
**********************************************************************************************************************************
********************************************************************************
**# Appendix Table 1: Home bias in vehicle market shares
********************************************************************************
use "data.dta", clear
tab origin
des brand*
bysort country year brand2: egen reg_cyb=total(reg), missing
bysort country year brand2: gen share_cyb=reg_cyb/reg_cy
keep if country=="France" | country=="Germany"
duplicates drop country year brand2, force
twoway scatter share_cyb brand2 if country=="France", by(year)
bysort country year: egen rank=rank(share_cyb)
br country year brand brand2 rank share_cyb  origin
//I select Citroen, Renault, Peugeot, VW,Audi and BMW
sum share_cyb if brand2==17 & country=="France"
sum share_cyb if brand2==71 & country=="France"
sum share_cyb if brand2==67 & country=="France"
sum share_cyb if brand2==90 & country=="France"
sum share_cyb if brand2==6 & country=="France"
sum share_cyb if brand2==7 & country=="France"
sum share_cyb if brand2==17 & country=="Germany"
sum share_cyb if brand2==71 & country=="Germany"
sum share_cyb if brand2==67 & country=="Germany"
sum share_cyb if brand2==90 & country=="Germany"
sum share_cyb if brand2==6 & country=="Germany"
sum share_cyb if brand2==7 & country=="Germany"


********************************************************************************
**# Appendix Table 2: Mean Quality by vehicle attributes
********************************************************************************
use "data.dta", clear
//by fuel type
tab fuelcat_harmonized2 [fw=reg], sum(wtp_jcy1)
//by number of engine cylinders
tab enginecyl
sum wtp_jcy1 [fw=reg] if enginecyl<=4
sum wtp_jcy1 [fw=reg] if enginecyl==5 | enginecyl==6
sum wtp_jcy1 [fw=reg] if enginecyl>=7 & enginecyl<.
//by body type
tab bodytype_harmonized4 [fw=reg], sum(wtp_jcy1)


********************************************************************************
**# Appendix Table 3: First-stage results for demand estimation
********************************************************************************
* Model 1:baseline regression
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(i.country2#i.year i.groupmodel2#i.bodytype_harmonized4    ///
	i.segment2#i.year) first ffirst cluster(group_mdltrim)
est store nl
matrix beta=e(b)'
matrix list beta

reghdfe price sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2    ///
	tax energycost1 loghw logw logsize ,   ///
	absorb( i.country2#i.year i.groupmodel2#i.bodytype_harmonized4 i.segment2#i.year) vce(cluster group_mdltrim)
est store stage1_price

reghdfe logshare_cyjinso sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2    ///
	tax energycost1 loghw logw logsize ,   ///
	absorb(i.country2#i.year i.groupmodel2#i.bodytype_harmonized4 i.segment2#i.year) vce(cluster group_mdltrim)
est store stage1_sharejinso

reghdfe logshare_cyoins sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2    ///
	tax energycost1 loghw logw logsize ,   ///
	absorb(i.country2#i.year i.groupmodel2#i.bodytype_harmonized4 i.segment2#i.year) vce(cluster group_mdltrim)
est store stage1_shareoins

esttab stage1_price stage1_sharejinso stage1_shareoins,b se stats(N F r2_a rmse) label  br compress
esttab stage1_price stage1_sharejinso stage1_shareoins,b(%7.3e) se(%7.3e) stats(N F r2_a rmse) nostar label compress br

est tab stage1_price stage1_sharejinso stage1_shareoins, stats(N F r2_a rmse)


********************************************************************************
**# Appendix Table 4: Cluster se by firm for DID regressions
********************************************************************************
use "data.dta", clear
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(i.country2#i.year i.groupmodel2#i.bodytype_harmonized4    ///
	 i.segment2#i.year) first ffirst cluster(group_mdltrim)
est store nl
matrix beta=e(b)'
matrix list beta
gen snl=e(sample)

//normalize the loghw and logw
drop loghw3 logw2
gen loghw3=loghw/abs(beta[1,1])
gen logw2=logw/abs(beta[1,1])
label var loghw3 "loghw (hp/kg) normalized by the absolute value of price coefficient"
label var logw2 "logw (ton) normalized by the absolute value of price coefficient"

* 1.0 main regressions using firm level stringnecy
foreach x in wtp_jcy1 loghw3 logw2 logprice {
	reghdfe `x' c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3   ///
		[fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year) vce(cluster group_mdltrim) 
	est store wf0_`x'
	count if e(sample)==1
	matrix n_wf0_`x'=r(N)
	test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
	matrix test_wf0_`x'=(r(drop),  r(df_r), r(F), r(df), r(p))'
}

* output
esttab wf0_wtp_jcy1 wf0_loghw3 wf0_logw2 wf0_logprice,     ///
	b(%7.3f) se(%7.3f) stats(N r2_a) br compress order(*.post#* stringency_* _cons) nostar nogap
estimates tab wf0_wtp_jcy1 wf0_loghw3 wf0_logw2, b(%7.3f) se(%7.3f) stats(N r2_a)

matrix A=(test_wf0_wtp_jcy1, test_wf0_loghw3, test_wf0_logw2, test_wf0_logprice)
matrix list A

matrix B=(n_wf0_wtp_jcy1, n_wf0_loghw3, n_wf0_logw2, n_wf0_logprice)
matrix list B

* Main DID regressions using firm level stringency, cluster at firm-year level 
tab company
unique company year //195
foreach x in wtp_jcy1 loghw3 logw2 logprice {
	reghdfe `x' c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3   ///
		[fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year) vce(cluster company year) 
	est store wf0_`x'
	count if e(sample)==1
	matrix n_wf0_`x'=r(N)
	test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
	matrix test_wf0_`x'=(r(drop),  r(df_r), r(F), r(df), r(p))'
}

* output
esttab wf0_wtp_jcy1 wf0_loghw3 wf0_logw2 wf0_logprice, b(%7.3f) se(%7.3f) stats(N r2_a) br compress order(*.post#* stringency_* ) star(* 0.1 ** 0.05 *** 0.01) nogap noline

matrix A=(test_wf0_wtp_jcy1, test_wf0_loghw3, test_wf0_logw2, test_wf0_logprice)
matrix list A

matrix B=(n_wf0_wtp_jcy1, n_wf0_loghw3, n_wf0_logw2, n_wf0_logprice)
matrix list B



********************************************************************************
**# Appendix Table 5-10: demand estimation, other specifications
********************************************************************************
* 2.1 NL_s: one-level nest: segment
sort country year vehicle
//original
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize    ///
	(price logshare_cyjins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff),   ///
	absorb(fe_cy_nls=i.country2#i.year fe_mb_nls=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_nls=i.segment2#i.year, resid(resid_nls)) first ffirst cluster(group_mdltrim)
est store nls
matrix beta_nls=e(b)'
matrix list beta_nls
gen gamma_jcy_nls=logshare_jcy-beta_nls[1,1]*price-beta_nls[2,1]*logshare_cyjins-beta_nls[3,1]*tax   ///
	-beta_nls[4,1]*energycost1-beta_nls[5,1]*loghw-beta_nls[6,1]*logw-fe_cy_nls
gen wtp_jcy_nls=gamma_jcy_nls/beta_nls[1,1]*-1
//add constant
gen dep_nls=logshare_jcy-beta_nls[1,1]*price-beta_nls[2,1]*logshare_cyjins
reghdfe dep_nls tax energycost1 loghw logw logsize  ,   ///
	absorb( fe_cy_nls2=i.country2#i.year fe_mb_nls2=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_nls2=i.segment2#i.year, resid(resid_nls2))  vce(cluster group_mdltrim)
est store nls2
matrix beta_nls2=e(b)'
matrix list beta_nls2
//gamma and wtp
sum fe_cy_nls
gen gamma_jcy_nls2=dep_nls-beta_nls2[1,1]*tax-beta_nls2[2,1]*energycost1-beta_nls2[3,1]*loghw   ///
	-beta_nls2[4,1]*logw-fe_cy_nls2
gen wtp_jcy_nls2=gamma_jcy_nls2/beta_nls[1,1]*-1
gen gamma3_jcy_nls2=dep_nls-beta_nls2[1,1]*tax-beta_nls2[2,1]*energycost1-beta_nls2[3,1]*loghw   ///
	-beta_nls2[4,1]*logw-fe_cy_nls2-beta_nls2[6,1]
gen wtp3_jcy_nls2=gamma3_jcy_nls2/beta_nls[1,1]*-1
	
bysort country year: egen gamma_cy_nls=mean(gamma_jcy_nls)
bysort year: egen gamma_y_nls=mean(gamma_jcy_nls)
bysort country year: egen gamma_cy_nls2=mean(gamma_jcy_nls2)
bysort year: egen gamma_y_nls2=mean(gamma_jcy_nls2)
bysort country year: egen gamma3_cy_nls2=mean(gamma3_jcy_nls2)
bysort year: egen gamma3_y_nls2=mean(gamma3_jcy_nls2)

bysort country year: egen wtp_cy_nls=mean(wtp_jcy_nls)
bysort year: egen wtp_y_nls=mean(wtp_jcy_nls)
bysort country year: egen wtp_cy_nls2=mean(wtp_jcy_nls2)
bysort year: egen wtp_y_nls2=mean(wtp_jcy_nls2)
bysort country year: egen wtp3_cy_nls2=mean(wtp3_jcy_nls2)
bysort year: egen wtp3_y_nls2=mean(wtp3_jcy_nls2)
	
//elasticity and wtp for attributes
gen elastic_price_nls=((1-share_cyjins)/(1-beta_nls[2,1])+(1-share_cys)*share_cyjins)*beta_nls[1,1]*price
sum elastic_price_nls,d
bysort country year: asgen avgela_price_nls=elastic_price_nls, weight(reg)
sum avgela_price_nls,d

matrix wtp_nls=beta_nls/abs(beta_nls[1,1])
matrix list wtp_nls	
	
	
	
	
* 2.2 NL_o: one-level nest: origin
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize    ///
	(price logshare_cyjino=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff   ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_nest_o sum_width_nest_o sum_height_nest_o    ///
	sum_enginecyl_nest_o),   ///
	absorb(fe_cy_nlo=i.country2#i.year fe_mb_nlo=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_nlo=i.segment2#i.year, resid(resid_nlo)) first ffirst cluster(group_mdltrim)
est store nlo	
matrix beta_nlo=e(b)'
matrix list beta_nlo

gen gamma_jcy_nlo=logshare_jcy-beta_nlo[1,1]*price-beta_nlo[2,1]*logshare_cyjino-beta_nlo[3,1]*tax   ///
	-beta_nlo[4,1]*energycost1-beta_nlo[5,1]*loghw-beta_nlo[6,1]*logw-fe_cy_nlo
gen wtp_jcy_nlo=gamma_jcy_nlo/beta_nlo[1,1]*-1

//add constant
gen dep_nlo=logshare_jcy-beta_nlo[1,1]*price-beta_nlo[2,1]*logshare_cyjino
reghdfe dep_nlo tax energycost1 loghw logw logsize,   ///
	absorb(fe_cy_nlo2=i.country2#i.year fe_mb_nlo2=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_nlo2=i.segment2#i.year, resid(resid_nlo2))  vce(cluster group_mdltrim)
est store nlo2	
matrix beta_nlo2=e(b)'
matrix list beta_nlo2

sum fe_cy_nlo
gen gamma_jcy_nlo2=dep_nlo-beta_nlo2[1,1]*tax-beta_nlo2[2,1]*energycost1-beta_nlo2[3,1]*loghw    ///
	-beta_nlo2[4,1]*logw-fe_cy_nlo2
gen gamma3_jcy_nlo2=dep_nlo-beta_nlo2[1,1]*tax-beta_nlo2[2,1]*energycost1-beta_nlo2[3,1]*loghw    ///
	-beta_nlo2[4,1]*logw-fe_cy_nlo2-beta_nlo2[6,1]
gen wtp_jcy_nlo2=gamma_jcy_nlo2/beta_nlo[1,1]*-1
gen wtp3_jcy_nlo2=gamma3_jcy_nlo2/beta_nlo[1,1]*-1

	
bysort country year: egen gamma_cy_nlo=mean(gamma_jcy_nlo)
bysort year: egen gamma_y_nlo=mean(gamma_jcy_nlo)
bysort country year: egen gamma_cy_nlo2=mean(gamma_jcy_nlo2)
bysort year: egen gamma_y_nlo2=mean(gamma_jcy_nlo2)
bysort country year: egen gamma3_cy_nlo2=mean(gamma3_jcy_nlo2)
bysort year: egen gamma3_y_nlo2=mean(gamma3_jcy_nlo2)

bysort country year: egen wtp_cy_nlo=mean(wtp_jcy_nlo)
bysort year: egen wtp_y_nlo=mean(wtp_jcy_nlo)
bysort country year: egen wtp_cy_nlo2=mean(wtp_jcy_nlo2)
bysort year: egen wtp_y_nlo2=mean(wtp_jcy_nlo2)
bysort country year: egen wtp3_cy_nlo2=mean(wtp3_jcy_nlo2)
bysort year: egen wtp3_y_nlo2=mean(wtp3_jcy_nlo2)

//elasticity and wtp for attributes
gen elastic_price_nlo=((1-share_cyjino)/(1-beta_nlo[2,1])+(1-share_cyo)*share_cyjino)*beta_nlo[1,1]*price
sum elastic_price_nlo,d
bysort country year: asgen avgela_price_nlo=elastic_price_nlo, weight(reg)
sum avgela_price_nlo,d

matrix wtp_nlo=beta_nlo/abs(beta_nlo[1,1])
matrix list wtp_nlo



* 2.3 NL_so:one-level nest: segment-origin
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize    ///
	(price logshare_cyjinso=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_nest_so sum_width_nest_so sum_height_nest_so    ///
	sum_enginecyl_nest_so),   ///
	absorb(fe_cy_nlso=i.country2#i.year fe_mb_nlso=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_nlso=i.segment2#i.year, resid(resid_nlso)) first ffirst cluster(group_mdltrim)
est store nlso	
matrix beta_nlso=e(b)'
matrix list beta_nlso

gen gamma_jcy_nlso=logshare_jcy-beta_nlso[1,1]*price-beta_nlso[2,1]*logshare_cyjinso-beta_nlso[3,1]*tax   ///
	-beta_nlso[4,1]*energycost1-beta_nlso[5,1]*loghw-beta_nlso[6,1]*logw-fe_cy_nlso	
gen wtp_jcy_nlso=gamma_jcy_nlso/beta_nlso[1,1]*-1

//add constant
gen dep_nlso=logshare_jcy-beta_nlso[1,1]*price-beta_nlso[2,1]*logshare_cyjinso
reghdfe dep_nlso tax energycost1 loghw logw logsize  ,   ///
	absorb( fe_cy_nlso2=i.country2#i.year fe_mb_nlso2=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_nlso2=i.segment2#i.year, resid(resid_nlso2))  vce(cluster group_mdltrim)
est store nlso2
matrix beta_nlso2=e(b)'
matrix list beta_nlso2

gen gamma_jcy_nlso2=dep_nlso-beta_nlso2[1,1]*tax-beta_nlso2[2,1]*energycost1-beta_nlso2[3,1]*loghw   ///
	-beta_nlso2[4,1]*logw-fe_cy_nlso2
gen wtp_jcy_nlso2=gamma_jcy_nlso2/beta_nlso[1,1]*-1
gen gamma3_jcy_nlso2=dep_nlso-beta_nlso2[1,1]*tax-beta_nlso2[2,1]*energycost1-beta_nlso2[3,1]*loghw   ///
	-beta_nlso2[4,1]*logw-fe_cy_nlso2-beta_nlso2[6,1]
gen wtp3_jcy_nlso2=gamma3_jcy_nlso2/beta_nlso[1,1]*-1
	
sum fe_cy_nlso

bysort country year: egen gamma_cy_nlso=mean(gamma_jcy_nlso)
bysort year: egen gamma_y_nlso=mean(gamma_jcy_nlso)
bysort country year: egen gamma_cy_nlso2=mean(gamma_jcy_nlso2)
bysort year: egen gamma_y_nlso2=mean(gamma_jcy_nlso2)
bysort country year: egen gamma3_cy_nlso2=mean(gamma3_jcy_nlso2)
bysort year: egen gamma3_y_nlso2=mean(gamma3_jcy_nlso2)

bysort country year: egen wtp_cy_nlso=mean(wtp_jcy_nlso)
bysort year: egen wtp_y_nlso=mean(wtp_jcy_nlso)
bysort country year: egen wtp_cy_nlso2=mean(wtp_jcy_nlso2)
bysort year: egen wtp_y_nlso2=mean(wtp_jcy_nlso2)
bysort country year: egen wtp3_cy_nlso2=mean(wtp3_jcy_nlso2)
bysort year: egen wtp3_y_nlso2=mean(wtp3_jcy_nlso2)

//elasticity and wtp for attributes
gen elastic_price_nlso=((1-share_cyjinso)/(1-beta_nlso[2,1])+(1-share_cyso)*share_cyjinso)*beta_nlso[1,1]*price
sum elastic_price_nlso,d
bysort country year: asgen avgela_price_nlso=elastic_price_nlso, weight(reg)
sum avgela_price_nlso,d

matrix wtp_nlso=beta_nlso/abs(beta_nlso[1,1])
matrix list wtp_nlso



* 2.4 model1 (main regression):two-level nest: segment-->origin, ivfe, model-bodytype FE, clustered mdltrim
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m1=i.country2#i.year fe_mb_m1=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_m1=i.segment2#i.year, resid(resid_m1)) first ffirst cluster(group_mdltrim)
est store m1
matrix beta_m1=e(b)'
matrix list beta_m1

gen gamma_jcy_m1=logshare_jcy-beta_m1[1,1]*price-beta_m1[2,1]*logshare_cyjinso-beta_m1[3,1]*logshare_cyoins    ///
	-beta_m1[4,1]*tax-beta_m1[5,1]*energycost1-beta_m1[6,1]*loghw-beta_m1[7,1]*logw-fe_cy_m1
gen wtp_jcy_m1=gamma_jcy_m1/beta_m1[1,1]*-1

//add constant
gen dep_m1=logshare_jcy-beta_m1[1,1]*price-beta_m1[2,1]*logshare_cyjinso-beta_m1[3,1]*logshare_cyoins  
reghdfe dep_m1 tax energycost1 loghw logw logsize  ,   ///
	absorb( fe_cy_m1t=i.country2#i.year fe_mb_m1t=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_m1t=i.segment2#i.year, resid(resid_m1t))  vce(cluster group_mdltrim)
est store m1t
matrix beta_m1t=e(b)'
matrix list beta_m1t

gen gamma_jcy_m1t=dep_m1-beta_m1t[1,1]*tax-beta_m1t[2,1]*energycost1-beta_m1t[3,1]*loghw-    ///
	beta_m1t[4,1]*logw-fe_cy_m1t
gen wtp_jcy_m1t=gamma_jcy_m1t/beta_m1[1,1]*-1
gen gamma3_jcy_m1t=dep_m1-beta_m1t[1,1]*tax-beta_m1t[2,1]*energycost1-beta_m1t[3,1]*loghw-    ///
	beta_m1t[4,1]*logw-fe_cy_m1t-beta_m1t[6,1]
gen wtp3_jcy_m1t=gamma3_jcy_m1t/beta_m1[1,1]*-1

sum fe_cy_m1

bysort country year: egen gamma_cy_m1=mean(gamma_jcy_m1)
bysort year: egen gamma_y_m1=mean(gamma_jcy_m1)
bysort country year: egen gamma_cy_m1t=mean(gamma_jcy_m1t)
bysort year: egen gamma_y_m1t=mean(gamma_jcy_m1t)
bysort country year: egen gamma3_cy_m1t=mean(gamma3_jcy_m1t)
bysort year: egen gamma3_y_m1t=mean(gamma3_jcy_m1t)

bysort country year: egen wtp_cy_m1=mean(wtp_jcy_m1)
bysort year: egen wtp_y_m1=mean(wtp_jcy_m1)
bysort country year: egen wtp_cy_m1t=mean(wtp_jcy_m1t)
bysort year: egen wtp_y_m1t=mean(wtp_jcy_m1t)
bysort country year: egen wtp3_cy_m1t=mean(wtp3_jcy_m1t)
bysort year: egen wtp3_y_m1t=mean(wtp3_jcy_m1t)

//elasticity and wtp for attributes
gen elastic_price_m1=((1-share_cyjinso)/(1-beta_m1[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m1[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m1[1,1]*price
sum elastic_price_m1,d
bysort country year: asgen avgela_price_m1=elastic_price_m1, weight(reg)
sum avgela_price_m1,d	

matrix wtp_m1=beta_m1/abs(beta_m1[1,1])
matrix list wtp_m1



* 2.5 model2:two-level nest: segment-->origin, ivfe, model-bodytype-trimlevel FE, clustered mdltrim
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize    ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m2=i.country2#i.year fe_mb_m2=i.groupmodel2#i.bodytype_harmonized4#i.trimlevel2    ///
	 fe_sy_m2=i.segment2#i.year, resid(resid_m2)) first ffirst cluster(group_mdltrim)
est store m2
matrix beta_m2=e(b)'
matrix list beta_m2

gen gamma_jcy_m2=logshare_jcy-beta_m2[1,1]*price-beta_m2[2,1]*logshare_cyjinso-beta_m2[3,1]*logshare_cyoins   ///
	-beta_m2[4,1]*tax-beta_m2[5,1]*energycost1-beta_m2[6,1]*loghw-beta_m2[7,1]*logw-fe_cy_m2
gen wtp_jcy_m2=gamma_jcy_m2/beta_m2[1,1]*-1

//add constant
sort country year vehicle
gen dep_m2=logshare_jcy-beta_m2[1,1]*price-beta_m2[2,1]*logshare_cyjinso-beta_m2[3,1]*logshare_cyoins
reghdfe dep_m2 tax energycost1 loghw logw logsize,    ///
	absorb(fe_cy_m2t=i.country2#i.year fe_mb_m2t=i.groupmodel2#i.bodytype_harmonized4#i.trimlevel2    ///
	 fe_sy_m2t=i.segment2#i.year, resid(resid_m2t)) vce(cluster group_mdltrim)
est store m2t
matrix beta_m2t=e(b)'
matrix list beta_m2t

gen gamma_jcy_m2t=dep_m2-beta_m2t[1,1]*tax-beta_m2t[2,1]*energycost1-beta_m2t[3,1]*loghw    ///
	-beta_m2t[4,1]*logw-fe_cy_m2t
gen wtp_jcy_m2t=gamma_jcy_m2t/beta_m2[1,1]*-1
gen gamma3_jcy_m2t=dep_m2-beta_m2t[1,1]*tax-beta_m2t[2,1]*energycost1-beta_m2t[3,1]*loghw    ///
	-beta_m2t[4,1]*logw-fe_cy_m2t-beta_m2t[6,1]
gen wtp3_jcy_m2t=gamma3_jcy_m2t/beta_m2[1,1]*-1

sum fe_cy_m2

bysort country year: egen gamma_cy_m2=mean(gamma_jcy_m2)
bysort year: egen gamma_y_m2=mean(gamma_jcy_m2)
bysort country year: egen gamma_cy_m2t=mean(gamma_jcy_m2t)
bysort year: egen gamma_y_m2t=mean(gamma_jcy_m2t)
bysort country year: egen gamma3_cy_m2t=mean(gamma3_jcy_m2t)
bysort year: egen gamma3_y_m2t=mean(gamma3_jcy_m2t)

bysort country year: egen wtp_cy_m2=mean(wtp_jcy_m2)
bysort year: egen wtp_y_m2=mean(wtp_jcy_m2)
bysort country year: egen wtp_cy_m2t=mean(wtp_jcy_m2t)
bysort year: egen wtp_y_m2t=mean(wtp_jcy_m2t)
bysort country year: egen wtp3_cy_m2t=mean(wtp3_jcy_m2t)
bysort year: egen wtp3_y_m2t=mean(wtp3_jcy_m2t)

//elasticity and wtp for attributes
gen elastic_price_m2=((1-share_cyjinso)/(1-beta_m2[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m2[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m2[1,1]*price
sum elastic_price_m2,d
bysort country year: asgen avgela_price_m2=elastic_price_m2, weight(reg)
sum avgela_price_m2,d	

matrix wtp_m2=beta_m2/abs(beta_m2[1,1])
matrix list wtp_m2





* 2.6 model4:two-level nest: segment-->origin, ivfe, model-bodytype FE, trimlevel FE, cluster(mdltrim)
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize    ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m4=i.country2#i.year fe_trim_m4=i.trimlevel2 fe_b_m4=i.bodytype_harmonized4    ///
	fe_mb_m4=i.groupmodel2#i.bodytype_harmonized4 fe_sy_m4=i.segment2#i.year, resid(resid_m4)) first ffirst cluster(group_mdltrim)	
est store m4
matrix beta_m4=e(b)'
matrix list beta_m4

gen gamma_jcy_m4=logshare_jcy-beta_m4[1,1]*price-beta_m4[2,1]*logshare_cyjinso-beta_m4[3,1]*logshare_cyoins    ///
	-beta_m4[4,1]*tax-beta_m4[5,1]*energycost1-beta_m4[6,1]*loghw-beta_m4[7,1]*logw-fe_cy_m4
gen wtp_jcy_m4=gamma_jcy_m4/beta_m4[1,1]*-1

//add constant
sort country year vehicle
gen dep_m4=logshare_jcy-beta_m4[1,1]*price-beta_m4[2,1]*logshare_cyjinso-beta_m4[3,1]*logshare_cyoins
reghdfe dep_m4 tax energycost1 loghw logw logsize,    ///
	absorb(fe_cy_m4t=i.country2#i.year fe_trim_m4t=i.trimlevel2 fe_b_m4t=i.bodytype_harmonized4    ///
	fe_mb_m4t=i.groupmodel2#i.bodytype_harmonized4 fe_sy_m4t=i.segment2#i.year, resid(resid_m4t)) vce(cluster group_mdltrim)	
est store m4t
matrix beta_m4t=e(b)'
matrix list beta_m4t

gen gamma_jcy_m4t=dep_m4-beta_m4t[1,1]*tax-beta_m4t[2,1]*energycost1-beta_m4t[3,1]*loghw    ///
	-beta_m4t[4,1]*logw-fe_cy_m4t
gen wtp_jcy_m4t=gamma_jcy_m4t/beta_m4[1,1]*-1
gen gamma3_jcy_m4t=dep_m4-beta_m4t[1,1]*tax-beta_m4t[2,1]*energycost1-beta_m4t[3,1]*loghw    ///
	-beta_m4t[4,1]*logw-fe_cy_m4t-beta_m4t[6,1]
gen wtp3_jcy_m4t=gamma3_jcy_m4t/beta_m4[1,1]*-1

sum fe_cy_m4

bysort country year: egen gamma_cy_m4=mean(gamma_jcy_m4)
bysort year: egen gamma_y_m4=mean(gamma_jcy_m4)
bysort country year: egen gamma_cy_m4t=mean(gamma_jcy_m4t)
bysort year: egen gamma_y_m4t=mean(gamma_jcy_m4t)
bysort country year: egen gamma3_cy_m4t=mean(gamma3_jcy_m4t)
bysort year: egen gamma3_y_m4t=mean(gamma3_jcy_m4t)

bysort country year: egen wtp_cy_m4=mean(wtp_jcy_m4)
bysort year: egen wtp_y_m4=mean(wtp_jcy_m4)
bysort country year: egen wtp_cy_m4t=mean(wtp_jcy_m4t)
bysort year: egen wtp_y_m4t=mean(wtp_jcy_m4t)
bysort country year: egen wtp3_cy_m4t=mean(wtp3_jcy_m4t)
bysort year: egen wtp3_y_m4t=mean(wtp3_jcy_m4t)


//elasticity and wtp for attributes
gen elastic_price_m4=((1-share_cyjinso)/(1-beta_m4[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m4[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m4[1,1]*price
sum elastic_price_m4,d
bysort country year: asgen avgela_price_m4=elastic_price_m4, weight(reg)
sum avgela_price_m4,d	

matrix wtp_m4=beta_m4/abs(beta_m4[1,1])
matrix list wtp_m4




* 2.7 model 5: nl_iv, nl for two nesting levels, iv, model-bodytype FE, cluster se, model-year FE
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m5=i.country2#i.year fe_mb_m5=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_my_m5=i.groupmodel2#i.year fe_sy_m5=i.segment2#i.year, resid(resid_m5))    ///
	first ffirst cluster(group_mdltrim)
est store m5
matrix beta_m5=e(b)'
matrix list beta_m5

gen gamma_jcy_m5=logshare_jcy-beta_m5[1,1]*price-beta_m5[2,1]*logshare_cyjinso-beta_m5[3,1]*logshare_cyoins    ///
	-beta_m5[4,1]*tax-beta_m5[5,1]*energycost1-beta_m5[6,1]*loghw-beta_m5[7,1]*logw-fe_cy_m5
gen wtp_jcy_m5=gamma_jcy_m5/beta_m5[1,1]*-1

//add constant
sort country year vehicle
gen dep_m5=logshare_jcy-beta_m5[1,1]*price-beta_m5[2,1]*logshare_cyjinso-beta_m5[3,1]*logshare_cyoins
reghdfe dep_m5 tax energycost1 loghw logw logsize,    ///
	absorb(fe_cy_m5t=i.country2#i.year fe_mb_m5t=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_my_m5t=i.groupmodel2#i.year fe_sy_m5t=i.segment2#i.year, resid(resid_m5t))  vce(cluster group_mdltrim)
est store m5t
matrix beta_m5t=e(b)'
matrix list beta_m5t

gen gamma_jcy_m5t=dep_m5-beta_m5t[1,1]*tax-beta_m5t[2,1]*energycost1-beta_m5t[3,1]*loghw    ///
	-beta_m5t[4,1]*logw-fe_cy_m5t
gen wtp_jcy_m5t=gamma_jcy_m5t/beta_m5[1,1]*-1
gen gamma3_jcy_m5t=dep_m5-beta_m5t[1,1]*tax-beta_m5t[2,1]*energycost1-beta_m5t[3,1]*loghw    ///
	-beta_m5t[4,1]*logw-fe_cy_m5t-beta_m5t[6,1]
gen wtp3_jcy_m5t=gamma3_jcy_m5t/beta_m5[1,1]*-1

sum fe_cy_m5

bysort country year: egen gamma_cy_m5=mean(gamma_jcy_m5)
bysort year: egen gamma_y_m5=mean(gamma_jcy_m5)
bysort country year: egen gamma_cy_m5t=mean(gamma_jcy_m5t)
bysort year: egen gamma_y_m5t=mean(gamma_jcy_m5t)
bysort country year: egen gamma3_cy_m5t=mean(gamma3_jcy_m5t)
bysort year: egen gamma3_y_m5t=mean(gamma3_jcy_m5t)

bysort country year: egen wtp_cy_m5=mean(wtp_jcy_m5)
bysort year: egen wtp_y_m5=mean(wtp_jcy_m5)
bysort country year: egen wtp_cy_m5t=mean(wtp_jcy_m5t)
bysort year: egen wtp_y_m5t=mean(wtp_jcy_m5t)
bysort country year: egen wtp3_cy_m5t=mean(wtp3_jcy_m5t)
bysort year: egen wtp3_y_m5t=mean(wtp3_jcy_m5t)

//elasticity and wtp for attributes
gen elastic_price_m5=((1-share_cyjinso)/(1-beta_m5[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m5[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m5[1,1]*price
sum elastic_price_m5,d
bysort country year: asgen avgela_price_m5=elastic_price_m5, weight(reg)
sum avgela_price_m5,d	

matrix wtp_m5=beta_m5/abs(beta_m5[1,1])
matrix list wtp_m5


* 2.8 model 6: nl_iv, model-bodytype FE, cluster se, congestion effect
* 1) calculate the number of products in each nest
bysort country year segment2 origin: gen num_so=_N
bysort country year segment2 origin: gen lognum_so=log(num_so)
bysort country year segment2: gen num_s=_N
bysort country year segment2: gen lognum_s=log(num_s)
sum num_so num_s lognum_so lognum_s

* 2) calculate the distance across K attributes between product j and all other products in the same nest
foreach x in price fuelcons_avg1 enginehorsepower grossvehicleweight size    ///
	bodytype_harmonized2 numberofdoors drivenwheels2 transtype fuelcat enginecyl numberofgears{
	bysort country year segment2 origin: egen mean_`x'_so=mean(`x')
	bysort country year segment2 origin: egen sd_`x'_so=sd(`x')
	bysort country year segment2: egen mean_`x'_s=mean(`x')
	bysort country year segment2: egen sd_`x'_s=sd(`x')
}
foreach x in price fuelcons_avg1 enginehorsepower grossvehicleweight size    ///
	bodytype_harmonized2 numberofdoors drivenwheels2 transtype fuelcat enginecyl numberofgears{
	gen dis_`x'_so=((`x'-mean_`x'_so)/sd_`x'_so)^2
	gen dis_`x'_s=((`x'-mean_`x'_s)/sd_`x'_s)^2
}	
egen dis_jinso=rowtotal(dis_*_so), missing
egen dis_jins=rowtotal(dis_*_s),missing
gen r_jinso=sqrt(dis_jinso)
gen r_jins=sqrt(dis_jins)
sum r_jinso r_jins,d    //I find r_jinso is similar to r_jins, this is because origin doesn't divide segment into many subgroups

* 3) model: two-level nest: segment-->origin, ivfe, model-bodytype FE, clustered mdltrim, num_so num_s, r_jinso r_jins
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize lognum_so r_jinso    ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m6=i.country2#i.year fe_mb_m6=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_m6=i.segment2#i.year, resid(resid_m6)) first ffirst cluster(group_mdltrim)
est store m6
matrix beta_m6=e(b)'
matrix list beta_m6

gen gamma_jcy_m6=logshare_jcy-beta_m6[1,1]*price-beta_m6[2,1]*logshare_cyjinso-beta_m6[3,1]*logshare_cyoins    ///
	-beta_m6[4,1]*tax-beta_m6[5,1]*energycost1-beta_m6[6,1]*loghw-beta_m6[7,1]*logw-fe_cy_m6    ///
	-beta_m6[9,1]*lognum_so-beta_m6[10,1]*r_jinso
gen wtp_jcy_m6=gamma_jcy_m6/beta_m6[1,1]*-1

//add constant
sort country year vehicle
gen dep_m6=logshare_jcy-beta_m6[1,1]*price-beta_m6[2,1]*logshare_cyjinso-beta_m6[3,1]*logshare_cyoins
reghdfe dep_m6 tax energycost1 loghw logw logsize lognum_so r_jinso,    ///
	absorb(fe_cy_m6t=i.country2#i.year fe_mb_m6t=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_m6t=i.segment2#i.year, resid(resid_m6t))  vce(cluster group_mdltrim)
est store m6t
matrix beta_m6t=e(b)'
matrix list beta_m6t

gen gamma_jcy_m6t=dep_m6-beta_m6t[1,1]*tax-beta_m6t[2,1]*energycost1-beta_m6t[3,1]*loghw    ///
	-beta_m6t[4,1]*logw-fe_cy_m6t-beta_m6t[6,1]*lognum_so-beta_m6t[7,1]*r_jinso
gen wtp_jcy_m6t=gamma_jcy_m6t/beta_m6[1,1]*-1
gen gamma3_jcy_m6t=dep_m6-beta_m6t[1,1]*tax-beta_m6t[2,1]*energycost1-beta_m6t[3,1]*loghw    ///
	-beta_m6t[4,1]*logw-fe_cy_m6t-beta_m6t[6,1]*lognum_so-beta_m6t[7,1]*r_jinso-beta_m6t[8,1]
gen wtp3_jcy_m6t=gamma3_jcy_m6t/beta_m6[1,1]*-1

sum fe_cy_m6

bysort country year: egen gamma_cy_m6=mean(gamma_jcy_m6)
bysort year: egen gamma_y_m6=mean(gamma_jcy_m6)
bysort country year: egen gamma_cy_m6t=mean(gamma_jcy_m6t)
bysort year: egen gamma_y_m6t=mean(gamma_jcy_m6t)
bysort country year: egen gamma3_cy_m6t=mean(gamma3_jcy_m6t)
bysort year: egen gamma3_y_m6t=mean(gamma3_jcy_m6t)

bysort country year: egen wtp_cy_m6=mean(wtp_jcy_m6)
bysort year: egen wtp_y_m6=mean(wtp_jcy_m6)
bysort country year: egen wtp_cy_m6t=mean(wtp_jcy_m6t)
bysort year: egen wtp_y_m6t=mean(wtp_jcy_m6t)
bysort country year: egen wtp3_cy_m6t=mean(wtp3_jcy_m6t)
bysort year: egen wtp3_y_m6t=mean(wtp3_jcy_m6t)

//elasticity and wtp for attributes
gen elastic_price_m6=((1-share_cyjinso)/(1-beta_m6[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m6[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m6[1,1]*price
sum elastic_price_m6,d
bysort country year: asgen avgela_price_m6=elastic_price_m6, weight(reg)
sum avgela_price_m6,d	

matrix wtp_m6=beta_m6/abs(beta_m6[1,1])
matrix list wtp_m6



* 2.9 model3: similar to model 1 but with robust se
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m3=i.country2#i.year fe_mb_m3=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_m3=i.segment2#i.year, resid(resid_m3)) first ffirst robust
est store m3
matrix beta_m3=e(b)'
matrix list beta_m3

gen gamma_jcy_m3=logshare_jcy-beta_m3[1,1]*price-beta_m3[2,1]*logshare_cyjinso-beta_m3[3,1]*logshare_cyoins    ///
	-beta_m3[4,1]*tax-beta_m3[5,1]*energycost1-beta_m3[6,1]*loghw-beta_m3[7,1]*logw-fe_cy_m3
gen wtp_jcy_m3=gamma_jcy_m3/beta_m3[1,1]*-1

//add constant
gen dep_m3=logshare_jcy-beta_m3[1,1]*price-beta_m3[2,1]*logshare_cyjinso-beta_m3[3,1]*logshare_cyoins  
reghdfe dep_m3 tax energycost1 loghw logw logsize  ,   ///
	absorb( fe_cy_m3t=i.country2#i.year fe_mb_m3t=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_m3t=i.segment2#i.year, resid(resid_m3t))  vce(robust)
est store m3t
matrix beta_m3t=e(b)'
matrix list beta_m3t

gen gamma_jcy_m3t=dep_m3-beta_m3t[1,1]*tax-beta_m3t[2,1]*energycost1-beta_m3t[3,1]*loghw-    ///
	beta_m3t[4,1]*logw-fe_cy_m3t
gen wtp_jcy_m3t=gamma_jcy_m3t/beta_m3[1,1]*-1
gen gamma3_jcy_m3t=dep_m3-beta_m3t[1,1]*tax-beta_m3t[2,1]*energycost1-beta_m3t[3,1]*loghw-    ///
	beta_m3t[4,1]*logw-fe_cy_m3t-beta_m3t[6,1]
gen wtp3_jcy_m3t=gamma3_jcy_m3t/beta_m3[1,1]*-1

sum fe_cy_m3

bysort country year: egen gamma_cy_m3=mean(gamma_jcy_m3)
bysort year: egen gamma_y_m3=mean(gamma_jcy_m3)
bysort country year: egen gamma_cy_m3t=mean(gamma_jcy_m3t)
bysort year: egen gamma_y_m3t=mean(gamma_jcy_m3t)
bysort country year: egen gamma3_cy_m3t=mean(gamma3_jcy_m3t)
bysort year: egen gamma3_y_m3t=mean(gamma3_jcy_m3t)

bysort country year: egen wtp_cy_m3=mean(wtp_jcy_m3)
bysort year: egen wtp_y_m3=mean(wtp_jcy_m3)
bysort country year: egen wtp_cy_m3t=mean(wtp_jcy_m3t)
bysort year: egen wtp_y_m3t=mean(wtp_jcy_m3t)
bysort country year: egen wtp3_cy_m3t=mean(wtp3_jcy_m3t)
bysort year: egen wtp3_y_m3t=mean(wtp3_jcy_m3t)


* 2.10 model 8: allow the preference parameters to vary across time
//original
ivreghdfe logshare_jcy tax i.year#c.energycost1 i.year#c.loghw i.year#c.logw i.year#c.logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m8=i.country2#i.year fe_mb_m8=i.groupmodel2#i.bodytype_harmonized4     ///
	 fe_sy_m8=i.segment2#i.year, resid(resid_m8)) first ffirst cluster(group_mdltrim)
est store m8
matrix beta_m8=e(b)'
matrix list beta_m8
esttab m8, b(%7.3f) se(%7.3f) stats(N r2_a) br compress nostar

gen gamma_jcy_m8=logshare_jcy-beta_m8[1,1]*price-beta_m8[2,1]*logshare_cyjinso-beta_m8[3,1]*logshare_cyoins -beta_m8[4,1]*tax        ///
	-beta_m8[5,1]*energycost1*cond(year==2005,1,0)-beta_m8[6,1]*energycost1*cond(year==2006,1,0)-beta_m8[7,1]*energycost1*cond(year==2007,1,0)-beta_m8[8,1]*energycost1*cond(year==2008,1,0)		///
	-beta_m8[9,1]*energycost1*cond(year==2009,1,0)-beta_m8[10,1]*energycost1*cond(year==2010,1,0)-beta_m8[11,1]*energycost1*cond(year==2011,1,0)-beta_m8[12,1]*energycost1*cond(year==2012,1,0)		///
	-beta_m8[13,1]*energycost1*cond(year==2013,1,0)-beta_m8[14,1]*energycost1*cond(year==2014,1,0)-beta_m8[15,1]*energycost1*cond(year==2015,1,0)-beta_m8[16,1]*energycost1*cond(year==2016,1,0)-beta_m8[17,1]*energycost1*cond(year==2017,1,0)		///
	-beta_m8[18,1]*loghw*cond(year==2005,1,0)-beta_m8[19,1]*loghw*cond(year==2006,1,0)-beta_m8[20,1]*loghw*cond(year==2007,1,0)-beta_m8[21,1]*loghw*cond(year==2008,1,0)		///
	-beta_m8[22,1]*loghw*cond(year==2009,1,0)-beta_m8[23,1]*loghw*cond(year==2010,1,0)-beta_m8[24,1]*loghw*cond(year==2011,1,0)-beta_m8[25,1]*loghw*cond(year==2012,1,0)		///
	-beta_m8[26,1]*loghw*cond(year==2013,1,0)-beta_m8[27,1]*loghw*cond(year==2014,1,0)-beta_m8[28,1]*loghw*cond(year==2015,1,0)-beta_m8[29,1]*loghw*cond(year==2016,1,0)-beta_m8[30,1]*loghw*cond(year==2017,1,0)		///
	-beta_m8[31,1]*logw*cond(year==2005,1,0)-beta_m8[32,1]*logw*cond(year==2006,1,0)-beta_m8[33,1]*logw*cond(year==2007,1,0)-beta_m8[34,1]*logw*cond(year==2008,1,0)		///
	-beta_m8[35,1]*logw*cond(year==2009,1,0)-beta_m8[36,1]*logw*cond(year==2010,1,0)-beta_m8[37,1]*logw*cond(year==2011,1,0)-beta_m8[38,1]*logw*cond(year==2012,1,0)		///
	-beta_m8[39,1]*logw*cond(year==2013,1,0)-beta_m8[40,1]*logw*cond(year==2014,1,0)-beta_m8[41,1]*logw*cond(year==2015,1,0)-beta_m8[42,1]*logw*cond(year==2016,1,0)-beta_m8[43,1]*logw*cond(year==2017,1,0)		///
	-fe_cy_m8
gen wtp_jcy_m8=gamma_jcy_m8/beta_m8[1,1]*-1

//add constant
gen dep_m8=logshare_jcy-beta_m8[1,1]*price-beta_m8[2,1]*logshare_cyjinso-beta_m8[3,1]*logshare_cyoins  
reghdfe dep_m8 tax i.year#c.energycost1 i.year#c.loghw i.year#c.logw i.year#c.logsize  ,   ///
	absorb( fe_cy_m8t=i.country2#i.year fe_mb_m8t=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_m8t=i.segment2#i.year, resid(resid_m8t))  vce(cluster group_mdltrim)
est store m8t
matrix beta_m8t=e(b)'
matrix list beta_m8t

gen gamma_jcy_m8t=dep_m8-beta_m8t[1,1]*tax		///
	-beta_m8t[2,1]*energycost1*cond(year==2005,1,0)-beta_m8t[3,1]*energycost1*cond(year==2006,1,0)-beta_m8t[4,1]*energycost1*cond(year==2007,1,0)-beta_m8t[5,1]*energycost1*cond(year==2008,1,0)		///
	-beta_m8t[6,1]*energycost1*cond(year==2009,1,0)-beta_m8t[7,1]*energycost1*cond(year==2010,1,0)-beta_m8t[8,1]*energycost1*cond(year==2011,1,0)-beta_m8t[9,1]*energycost1*cond(year==2012,1,0)		///
	-beta_m8t[10,1]*energycost1*cond(year==2013,1,0)-beta_m8t[11,1]*energycost1*cond(year==2014,1,0)-beta_m8t[12,1]*energycost1*cond(year==2015,1,0)-beta_m8t[13,1]*energycost1*cond(year==2016,1,0)-beta_m8t[14,1]*energycost1*cond(year==2017,1,0)		///
	-beta_m8t[15,1]*loghw*cond(year==2005,1,0)-beta_m8t[16,1]*loghw*cond(year==2006,1,0)-beta_m8t[17,1]*loghw*cond(year==2007,1,0)-beta_m8t[18,1]*loghw*cond(year==2008,1,0)		///
	-beta_m8t[19,1]*loghw*cond(year==2009,1,0)-beta_m8t[20,1]*loghw*cond(year==2010,1,0)-beta_m8t[21,1]*loghw*cond(year==2011,1,0)-beta_m8t[22,1]*loghw*cond(year==2012,1,0)		///
	-beta_m8t[23,1]*loghw*cond(year==2013,1,0)-beta_m8t[24,1]*loghw*cond(year==2014,1,0)-beta_m8t[25,1]*loghw*cond(year==2015,1,0)-beta_m8t[26,1]*loghw*cond(year==2016,1,0)-beta_m8t[27,1]*loghw*cond(year==2017,1,0)		///
	-beta_m8t[28,1]*logw*cond(year==2005,1,0)-beta_m8t[29,1]*logw*cond(year==2006,1,0)-beta_m8t[30,1]*logw*cond(year==2007,1,0)-beta_m8t[31,1]*logw*cond(year==2008,1,0)		///
	-beta_m8t[32,1]*logw*cond(year==2009,1,0)-beta_m8t[33,1]*logw*cond(year==2010,1,0)-beta_m8t[34,1]*logw*cond(year==2011,1,0)-beta_m8t[35,1]*logw*cond(year==2012,1,0)		///
	-beta_m8t[36,1]*logw*cond(year==2013,1,0)-beta_m8t[37,1]*logw*cond(year==2014,1,0)-beta_m8t[38,1]*logw*cond(year==2015,1,0)-beta_m8t[39,1]*logw*cond(year==2016,1,0)-beta_m8t[40,1]*logw*cond(year==2017,1,0)		///
	-fe_cy_m8t
gen wtp_jcy_m8t=gamma_jcy_m8t/beta_m8[1,1]*-1
gen gamma3_jcy_m8t=dep_m8-beta_m8t[1,1]*tax		///
	-beta_m8t[2,1]*energycost1*cond(year==2005,1,0)-beta_m8t[3,1]*energycost1*cond(year==2006,1,0)-beta_m8t[4,1]*energycost1*cond(year==2007,1,0)-beta_m8t[5,1]*energycost1*cond(year==2008,1,0)		///
	-beta_m8t[6,1]*energycost1*cond(year==2009,1,0)-beta_m8t[7,1]*energycost1*cond(year==2010,1,0)-beta_m8t[8,1]*energycost1*cond(year==2011,1,0)-beta_m8t[9,1]*energycost1*cond(year==2012,1,0)		///
	-beta_m8t[10,1]*energycost1*cond(year==2013,1,0)-beta_m8t[11,1]*energycost1*cond(year==2014,1,0)-beta_m8t[12,1]*energycost1*cond(year==2015,1,0)-beta_m8t[13,1]*energycost1*cond(year==2016,1,0)-beta_m8t[14,1]*energycost1*cond(year==2017,1,0)		///
	-beta_m8t[15,1]*loghw*cond(year==2005,1,0)-beta_m8t[16,1]*loghw*cond(year==2006,1,0)-beta_m8t[17,1]*loghw*cond(year==2007,1,0)-beta_m8t[18,1]*loghw*cond(year==2008,1,0)		///
	-beta_m8t[19,1]*loghw*cond(year==2009,1,0)-beta_m8t[20,1]*loghw*cond(year==2010,1,0)-beta_m8t[21,1]*loghw*cond(year==2011,1,0)-beta_m8t[22,1]*loghw*cond(year==2012,1,0)		///
	-beta_m8t[23,1]*loghw*cond(year==2013,1,0)-beta_m8t[24,1]*loghw*cond(year==2014,1,0)-beta_m8t[25,1]*loghw*cond(year==2015,1,0)-beta_m8t[26,1]*loghw*cond(year==2016,1,0)-beta_m8t[27,1]*loghw*cond(year==2017,1,0)		///
	-beta_m8t[28,1]*logw*cond(year==2005,1,0)-beta_m8t[29,1]*logw*cond(year==2006,1,0)-beta_m8t[30,1]*logw*cond(year==2007,1,0)-beta_m8t[31,1]*logw*cond(year==2008,1,0)		///
	-beta_m8t[32,1]*logw*cond(year==2009,1,0)-beta_m8t[33,1]*logw*cond(year==2010,1,0)-beta_m8t[34,1]*logw*cond(year==2011,1,0)-beta_m8t[35,1]*logw*cond(year==2012,1,0)		///
	-beta_m8t[36,1]*logw*cond(year==2013,1,0)-beta_m8t[37,1]*logw*cond(year==2014,1,0)-beta_m8t[38,1]*logw*cond(year==2015,1,0)-beta_m8t[39,1]*logw*cond(year==2016,1,0)-beta_m8t[40,1]*logw*cond(year==2017,1,0)		///
	-fe_cy_m8t-beta_m8t[54,1]
gen wtp3_jcy_m8t=gamma3_jcy_m8t/beta_m8[1,1]*-1

bysort country year: egen gamma_cy_m8=mean(gamma_jcy_m8)
bysort year: egen gamma_y_m8=mean(gamma_jcy_m8)
bysort country year: egen gamma_cy_m8t=mean(gamma_jcy_m8t)
bysort year: egen gamma_y_m8t=mean(gamma_jcy_m8t)
bysort country year: egen gamma3_cy_m8t=mean(gamma3_jcy_m8t)
bysort year: egen gamma3_y_m8t=mean(gamma3_jcy_m8t)

bysort country year: egen wtp_cy_m8=mean(wtp_jcy_m8)
bysort year: egen wtp_y_m8=mean(wtp_jcy_m8)
bysort country year: egen wtp_cy_m8t=mean(wtp_jcy_m8t)
bysort year: egen wtp_y_m8t=mean(wtp_jcy_m8t)
bysort country year: egen wtp3_cy_m8t=mean(wtp3_jcy_m8t)
bysort year: egen wtp3_y_m8t=mean(wtp3_jcy_m8t)

//elasticity and wtp for attributes
gen elastic_price_m8=((1-share_cyjinso)/(1-beta_m8[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m8[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m8[1,1]*price
sum elastic_price_m8,d
bysort country year: asgen avgela_price_m8=elastic_price_m8, weight(reg)
sum avgela_price_m8,d	

matrix wtp_m8=beta_m8/abs(beta_m8[1,1])
matrix list wtp_m8





* 2.11 model 9: other demand estimation, Dummy for luxury cars
use "data2.dta", clear
* 1.generate dummy variable for luxury cars
sum price [fw=reg], d
gen routlier3_price=cond(price>r(p99),1,0)
sum price [fw=reg], d
gen routlier4_price=cond(price>r(p95),1,0)

sum enginehorsepower grossvehicleweight loghw loghw2
drop loghw
gen loghw=log(enginehorsepower*100/(grossvehicleweight*1000))
label var loghw "log(enginehorsepower*100/(grossvehicleweight*1000))"

count
sum price

ivreghdfe logshare_jcy tax energycost1 loghw logw logsize routlier3_price  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_d1=i.country2#i.year fe_mb_d1=i.groupmodel2#i.bodytype_harmonized4 		///
	 fe_sy_d1=i.segment2#i.year, resid(resid_d1)) first ffirst cluster(group_mdltrim)
est store nl_d1
matrix beta_d1=e(b)'
matrix list beta_d1   

predict fe_nl_d1, d

//add constant
matrix list beta_d1
gen dep_nl_d1=logshare_jcy-beta_d1[1,1]*price-beta_d1[2,1]*logshare_cyjinso-beta_d1[3,1]*logshare_cyoins
reghdfe dep_nl_d1 tax energycost1 loghw logw logsize routlier3_price  ,   ///
	absorb(fe_cy_d1t=i.country2#i.year fe_mb_d1t=i.groupmodel2#i.bodytype_harmonized4 		///
	 fe_sy_d1t=i.segment2#i.year, resid(resid_d1t)) vce(cluster group_mdltrim)
est store nl_d1t
matrix beta_d1t=e(b)'
matrix list beta_d1t
predict fe_nl_d1t, d
sum fe_cy_d1   

//gamma and wtp
gen gamma_jcy_d1=logshare_jcy-beta_d1[1,1]*price-beta_d1[2,1]*logshare_cyjinso-beta_d1[3,1]*logshare_cyoins    ///
	-beta_d1[4,1]*tax-beta_d1[5,1]*energycost1-beta_d1[6,1]*loghw-beta_d1[7,1]*logw-fe_cy_d1-beta_d1[9,1]*routlier3_price
gen wtp_jcy_d1=gamma_jcy_d1/beta_d1[1,1]*-1

gen gamma_jcy_d1t=dep_nl_d1-beta_d1t[1,1]*tax-beta_d1t[2,1]*energycost1-beta_d1t[3,1]*loghw-    ///
	beta_d1t[4,1]*logw-fe_cy_d1t-beta_d1t[6,1]*routlier3_price
gen gamma3_jcy_d1t=dep_nl_d1-beta_d1t[1,1]*tax-beta_d1t[2,1]*energycost1-beta_d1t[3,1]*loghw-    ///
	beta_d1t[4,1]*logw-fe_cy_d1t-beta_d1t[6,1]*routlier3_price-beta_d1t[7,1]
gen wtp_jcy_d1t=gamma_jcy_d1t/beta_d1[1,1]*-1
gen wtp3_jcy_d1t=gamma3_jcy_d1t/beta_d1[1,1]*-1


//elasticity and wtp for attributes
gen elastic_price_d1=((1-share_cyjinso)/(1-beta_d1[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_d1[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_d1[1,1]*price
sum elastic_price_d1,d
bysort country year: asgen avgela_price_d1=elastic_price_d1, weight(reg)
sum avgela_price_d1,d	

matrix wtp_d1=beta_d1/abs(beta_d1[1,1])
matrix list wtp_d1

//quality index by country-year
bysort country year: egen gamma_cy_d1=mean(gamma_jcy_d1)
bysort country year: egen wtp_cy_d1=mean(wtp_jcy_d1)
bysort country year: egen gamma_cy_d1t=mean(gamma_jcy_d1t)
bysort country year: egen wtp_cy_d1t=mean(wtp_jcy_d1t)
bysort year: egen wtp_y_d1=mean(wtp_jcy_d1)
bysort year: egen wtp_y_d1t=mean(wtp_jcy_d1t)

bysort country year: egen gamma3_cy_d1=mean(gamma3_jcy_d1)
bysort country year: egen wtp3_cy_d1=mean(wtp3_jcy_d1)
bysort country year: egen gamma3_cy_d1t=mean(gamma3_jcy_d1t)
bysort country year: egen wtp3_cy_d1t=mean(wtp3_jcy_d1t)
bysort year: egen wtp3_y_d1=mean(wtp3_jcy_d1)
bysort year: egen wtp3_y_d1t=mean(wtp3_jcy_d1t)

save "demand_luxury_20200830.dta", replace


*************************************
* output for appendix table 
*************************************
//appendix table 5: other demand structure
esttab m1 nls nlo nlso, b(%7.3f) se(%7.3f) stats(N r2_a) br compress nostar   ///
	order(price logshare_cyjins logshare_cyjino logshare_cyjinso logshare_cyoins tax energycost1 loghw logw logsize _cons)

//appendix table 8: other identification
esttab m1 m2 m4 m5 m6 m3, b(%7.3f) se(%7.3f) stats(N r2_a) br compress nostar   ///
	order(price logshare_cyjinso logshare_cyoins lognum_so r_jinso tax energycost1 loghw logw logsize)

//appendix table 8: varying parameters, identification model 8
esttab m8, b(%7.3f) se(%7.3f) stats(N r2_a) br compress nostar   ///
	order(price logshare_cyjinso logshare_cyoins tax energycost1 loghw logw logsize)

//appendix table 8: a dummy variable for luxury cars
esttab nl_d1, b(%7.3f) se(%7.3f) stats(N r2_a) br compress nostar


*************************************
* output for average price elasticity 
*************************************
//appendix table 6
sum avgela_price_m1, d 
matrix A=(r(p50)\r(mean)\r(sd)\r(p5)\r(p95))
foreach x in nls nlo nlso{
    sum avgela_price_`x', d
	matrix A=(A,(r(p50)\r(mean)\r(sd)\r(p5)\r(p95)))
}
matrix list A
matrix drop A
//appendix table 9
sum avgela_price_m1, d 
matrix A=(r(p50)\r(mean)\r(sd)\r(p5)\r(p95))
foreach x in m2 m4 m5 m6 d1 m8 m1{
    sum avgela_price_`x', d
	matrix A=(A,(r(p50)\r(mean)\r(sd)\r(p5)\r(p95)))
}
matrix list A
matrix drop A




********************************************************************************
**# Appendix Table 11-13 : demand estimation, different fixed effects
********************************************************************************
use "data.dta", clear
* 1.1 model1 (main regression):two-level nest: segment-->origin, ivfe, model-bodytype FE, clustered mdltrim
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m1=i.country2#i.year fe_mb_m1=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_m1=i.segment2#i.year, resid(resid_m1)) first ffirst cluster(group_mdltrim)
est store m1
matrix beta_m1=e(b)'
matrix list beta_m1

gen gamma_jcy_m1=logshare_jcy-beta_m1[1,1]*price-beta_m1[2,1]*logshare_cyjinso-beta_m1[3,1]*logshare_cyoins    ///
	-beta_m1[4,1]*tax-beta_m1[5,1]*energycost1-beta_m1[6,1]*loghw-beta_m1[7,1]*logw-fe_cy_m1
gen wtp_jcy_m1=gamma_jcy_m1/beta_m1[1,1]*-1

//elasticity and wtp for attributes
gen elastic_price_m1=((1-share_cyjinso)/(1-beta_m1[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m1[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m1[1,1]*price
sum elastic_price_m1,d
bysort country year: asgen avgela_price_m1=elastic_price_m1, weight(reg)
sum avgela_price_m1,d	



* 1.2 model2:two-level nest: segment-->origin, ivfe, model-bodytype-trimlevel FE, segment-year
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize    ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m2=i.country2#i.year fe_mb_m2=i.groupmodel2#i.bodytype_harmonized4#i.trimlevel2    ///
	 fe_sy_m2=i.segment2#i.year, resid(resid_m2)) first ffirst cluster(group_mdltrim)
est store m2
matrix beta_m2=e(b)'
matrix list beta_m2

gen gamma_jcy_m2=logshare_jcy-beta_m2[1,1]*price-beta_m2[2,1]*logshare_cyjinso-beta_m2[3,1]*logshare_cyoins   ///
	-beta_m2[4,1]*tax-beta_m2[5,1]*energycost1-beta_m2[6,1]*loghw-beta_m2[7,1]*logw-fe_cy_m2
gen wtp_jcy_m2=gamma_jcy_m2/beta_m2[1,1]*-1

//elasticity and wtp for attributes
gen elastic_price_m2=((1-share_cyjinso)/(1-beta_m2[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m2[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m2[1,1]*price
sum elastic_price_m2,d
bysort country year: asgen avgela_price_m2=elastic_price_m2, weight(reg)
sum avgela_price_m2,d	


* 1.3 model 3: two-level nest: segment-->origin, ivfe, model-bodytype-trimlevel-segment FE, segment-year
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m3=i.country2#i.year fe_mb_m3=i.groupmodel2#i.bodytype_harmonized4#i.trimlevel2#i.segment2   ///
	 fe_sy_m3=i.segment2#i.year, resid(resid_m3)) first ffirst cluster(group_mdltrim)
est store m3
matrix beta_m3=e(b)'
matrix list beta_m3

gen gamma_jcy_m3=logshare_jcy-beta_m3[1,1]*price-beta_m3[2,1]*logshare_cyjinso-beta_m3[3,1]*logshare_cyoins    ///
	-beta_m3[4,1]*tax-beta_m3[5,1]*energycost1-beta_m3[6,1]*loghw-beta_m3[7,1]*logw-fe_cy_m3
gen wtp_jcy_m3=gamma_jcy_m3/beta_m3[1,1]*-1

//elasticity and wtp for attributes
gen elastic_price_m3=((1-share_cyjinso)/(1-beta_m3[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m3[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m3[1,1]*price
sum elastic_price_m3,d
bysort country year: asgen avgela_price_m3=elastic_price_m3, weight(reg)
sum avgela_price_m3,d	


* 1.4 model 4: two-level nest: segment-->origin, ivfe, model-bodytype-trimlevel-segment-fueltype FE, segment-year
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m4=i.country2#i.year fe_mb_m4=i.groupmodel2#i.bodytype_harmonized4#i.trimlevel2#i.segment2#i.fuelcat_harmonized2   ///
	 fe_sy_m4=i.segment2#i.year, resid(resid_m4)) first ffirst cluster(group_mdltrim)
est store m4
matrix beta_m4=e(b)'
matrix list beta_m4

gen gamma_jcy_m4=logshare_jcy-beta_m4[1,1]*price-beta_m4[2,1]*logshare_cyjinso-beta_m4[3,1]*logshare_cyoins    ///
	-beta_m4[4,1]*tax-beta_m4[5,1]*energycost1-beta_m4[6,1]*loghw-beta_m4[7,1]*logw-fe_cy_m4
gen wtp_jcy_m4=gamma_jcy_m4/beta_m4[1,1]*-1

//elasticity and wtp for attributes
gen elastic_price_m4=((1-share_cyjinso)/(1-beta_m4[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m4[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m4[1,1]*price
sum elastic_price_m4,d
bysort country year: asgen avgela_price_m4=elastic_price_m4, weight(reg)
sum avgela_price_m4,d	

* 1.5 model 5: two-level nest: segment-->origin, ivfe, vehicle FE, segment-year
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m5=i.country2#i.year fe_mb_m5=vehicle   ///
	 fe_sy_m5=i.segment2#i.year, resid(resid_m5)) first ffirst cluster(group_mdltrim)
est store m5
matrix beta_m5=e(b)'
matrix list beta_m5

gen gamma_jcy_m5=logshare_jcy-beta_m5[1,1]*price-beta_m5[2,1]*logshare_cyjinso-beta_m5[3,1]*logshare_cyoins    ///
	-beta_m5[4,1]*tax-beta_m5[5,1]*energycost1-beta_m5[6,1]*loghw-beta_m5[7,1]*logw-fe_cy_m5
gen wtp_jcy_m5=gamma_jcy_m5/beta_m5[1,1]*-1

//elasticity and wtp for attributes
gen elastic_price_m5=((1-share_cyjinso)/(1-beta_m5[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m5[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m5[1,1]*price
sum elastic_price_m5,d
bysort country year: asgen avgela_price_m5=elastic_price_m5, weight(reg)
sum avgela_price_m5,d	

* 1.6 model 6:model-bodytype FE, segment-year, groupmodel2-year, cy
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m6=i.country2#i.year fe_mb_m6=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_m6=i.segment2#i.year fe_my6=i.groupmodel2#i.year, resid(resid_m6)) first ffirst cluster(group_mdltrim)
est store m6
matrix beta_m6=e(b)'
matrix list beta_m6

gen gamma_jcy_m6=logshare_jcy-beta_m6[1,1]*price-beta_m6[2,1]*logshare_cyjinso-beta_m6[3,1]*logshare_cyoins    ///
	-beta_m6[4,1]*tax-beta_m6[5,1]*energycost1-beta_m6[6,1]*loghw-beta_m6[7,1]*logw-fe_cy_m6
gen wtp_jcy_m6=gamma_jcy_m6/beta_m6[1,1]*-1

//elasticity and wtp for attributes
gen elastic_price_m6=((1-share_cyjinso)/(1-beta_m6[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m6[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m6[1,1]*price
sum elastic_price_m6,d
bysort country year: asgen avgela_price_m6=elastic_price_m6, weight(reg)
sum avgela_price_m6,d	

* 1.7 model 7:model-bodytype FE, segment-year, bodytype-year, cy
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m7=i.country2#i.year fe_mb_m7=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_m7=i.segment2#i.year fe_by7=i.bodytype_harmonized4#i.year, resid(resid_m7)) first ffirst cluster(group_mdltrim)
est store m7
matrix beta_m7=e(b)'
matrix list beta_m7

gen gamma_jcy_m7=logshare_jcy-beta_m7[1,1]*price-beta_m7[2,1]*logshare_cyjinso-beta_m7[3,1]*logshare_cyoins    ///
	-beta_m7[4,1]*tax-beta_m7[5,1]*energycost1-beta_m7[6,1]*loghw-beta_m7[7,1]*logw-fe_cy_m7
gen wtp_jcy_m7=gamma_jcy_m7/beta_m7[1,1]*-1

//elasticity and wtp for attributes
gen elastic_price_m7=((1-share_cyjinso)/(1-beta_m7[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m7[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m7[1,1]*price
sum elastic_price_m7,d
bysort country year: asgen avgela_price_m7=elastic_price_m7, weight(reg)
sum avgela_price_m7,d	

* 1.8 model 8:model-bodytype FE, segment-year, trimlevel-year, cy
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m8=i.country2#i.year fe_mb_m8=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_m8=i.segment2#i.year fe_my8=i.trimlevel2#i.year, resid(resid_m8)) first ffirst cluster(group_mdltrim)
est store m8
matrix beta_m8=e(b)'
matrix list beta_m8

gen gamma_jcy_m8=logshare_jcy-beta_m8[1,1]*price-beta_m8[2,1]*logshare_cyjinso-beta_m8[3,1]*logshare_cyoins    ///
	-beta_m8[4,1]*tax-beta_m8[5,1]*energycost1-beta_m8[6,1]*loghw-beta_m8[7,1]*logw-fe_cy_m8
gen wtp_jcy_m8=gamma_jcy_m8/beta_m8[1,1]*-1

//elasticity and wtp for attributes
gen elastic_price_m8=((1-share_cyjinso)/(1-beta_m8[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m8[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m8[1,1]*price
sum elastic_price_m8,d
bysort country year: asgen avgela_price_m8=elastic_price_m8, weight(reg)
sum avgela_price_m8,d	

* 1.9 model 9:model-bodytype FE, segment-year, fuelcat-year, cy
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m9=i.country2#i.year fe_mb_m9=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_m9=i.segment2#i.year fe_fy9=i.fuelcat_harmonized2#i.year, resid(resid_m9)) first ffirst cluster(group_mdltrim)
est store m9
matrix beta_m9=e(b)'
matrix list beta_m9

gen gamma_jcy_m9=logshare_jcy-beta_m9[1,1]*price-beta_m9[2,1]*logshare_cyjinso-beta_m9[3,1]*logshare_cyoins    ///
	-beta_m9[4,1]*tax-beta_m9[5,1]*energycost1-beta_m9[6,1]*loghw-beta_m9[7,1]*logw-fe_cy_m9
gen wtp_jcy_m9=gamma_jcy_m9/beta_m9[1,1]*-1

//elasticity and wtp for attributes
gen elastic_price_m9=((1-share_cyjinso)/(1-beta_m9[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m9[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m9[1,1]*price
sum elastic_price_m9,d
bysort country year: asgen avgela_price_m9=elastic_price_m9, weight(reg)
sum avgela_price_m9,d	


* 1.10 model 10:model-bodytype FE, segment-year, groupmodel2-bodytype-year, cy
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m10=i.country2#i.year fe_mb_m10=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_m10=i.segment2#i.year fe_mby10=i.groupmodel2#i.bodytype_harmonized4#i.year, resid(resid_m10)) first ffirst cluster(group_mdltrim)
est store m10
matrix beta_m10=e(b)'
matrix list beta_m10

gen gamma_jcy_m10=logshare_jcy-beta_m10[1,1]*price-beta_m10[2,1]*logshare_cyjinso-beta_m10[3,1]*logshare_cyoins    ///
	-beta_m10[4,1]*tax-beta_m10[5,1]*energycost1-beta_m10[6,1]*loghw-beta_m10[7,1]*logw-fe_cy_m10
gen wtp_jcy_m10=gamma_jcy_m10/beta_m10[1,1]*-1

//elasticity and wtp for attributes
gen elastic_price_m10=((1-share_cyjinso)/(1-beta_m10[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m10[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m10[1,1]*price
sum elastic_price_m10,d
bysort country year: asgen avgela_price_m10=elastic_price_m10, weight(reg)
sum avgela_price_m10,d	

* 1.11 model 11:model-bodytype FE, segment-year, groupmodel2-year, bodytype-year, cy
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m11=i.country2#i.year fe_mb_m11=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_m11=i.segment2#i.year fe_my11=i.groupmodel2#i.year fe_by11=i.bodytype_harmonized4#i.year, resid(resid_m11)) first ffirst cluster(group_mdltrim)
est store m11
matrix beta_m11=e(b)'
matrix list beta_m11

gen gamma_jcy_m11=logshare_jcy-beta_m11[1,1]*price-beta_m11[2,1]*logshare_cyjinso-beta_m11[3,1]*logshare_cyoins    ///
	-beta_m11[4,1]*tax-beta_m11[5,1]*energycost1-beta_m11[6,1]*loghw-beta_m11[7,1]*logw-fe_cy_m11
gen wtp_jcy_m11=gamma_jcy_m11/beta_m11[1,1]*-1

//elasticity and wtp for attributes
gen elastic_price_m11=((1-share_cyjinso)/(1-beta_m11[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m11[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m11[1,1]*price
sum elastic_price_m11,d
bysort country year: asgen avgela_price_m11=elastic_price_m11, weight(reg)
sum avgela_price_m11,d	

* 1.12 model 12:model-bodytype FE, segment-year, groupmodel2-year, bodytype-year, trimlevel-year, cy
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m12=i.country2#i.year fe_mb_m12=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_m12=i.segment2#i.year fe_my12=i.groupmodel2#i.year fe_by12=i.bodytype_harmonized4#i.year fe_ty12=i.trimlevel2#i.year, resid(resid_m12)) first ffirst cluster(group_mdltrim)
est store m12
matrix beta_m12=e(b)'
matrix list beta_m12

gen gamma_jcy_m12=logshare_jcy-beta_m12[1,1]*price-beta_m12[2,1]*logshare_cyjinso-beta_m12[3,1]*logshare_cyoins    ///
	-beta_m12[4,1]*tax-beta_m12[5,1]*energycost1-beta_m12[6,1]*loghw-beta_m12[7,1]*logw-fe_cy_m12
gen wtp_jcy_m12=gamma_jcy_m12/beta_m12[1,1]*-1

//elasticity and wtp for attributes
gen elastic_price_m12=((1-share_cyjinso)/(1-beta_m12[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m12[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m12[1,1]*price
sum elastic_price_m12,d
bysort country year: asgen avgela_price_m12=elastic_price_m12, weight(reg)
sum avgela_price_m12,d	

* 1.13 model 13:model-bodytype FE, segment-year, groupmodel2-year, bodytype-year, trimlevel-year, fuelcat-year, cy
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m13=i.country2#i.year fe_mb_m13=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_m13=i.segment2#i.year fe_my13=i.groupmodel2#i.year fe_by13=i.bodytype_harmonized4#i.year fe_ty13=i.trimlevel2#i.year fe_fy13=i.fuelcat_harmonized2#i.year, resid(resid_m13)) first ffirst cluster(group_mdltrim)
est store m13
matrix beta_m13=e(b)'
matrix list beta_m13

gen gamma_jcy_m13=logshare_jcy-beta_m13[1,1]*price-beta_m13[2,1]*logshare_cyjinso-beta_m13[3,1]*logshare_cyoins    ///
	-beta_m13[4,1]*tax-beta_m13[5,1]*energycost1-beta_m13[6,1]*loghw-beta_m13[7,1]*logw-fe_cy_m13
gen wtp_jcy_m13=gamma_jcy_m13/beta_m13[1,1]*-1

//elasticity and wtp for attributes
gen elastic_price_m13=((1-share_cyjinso)/(1-beta_m13[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m13[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m13[1,1]*price
sum elastic_price_m13,d
bysort country year: asgen avgela_price_m13=elastic_price_m13, weight(reg)
sum avgela_price_m13,d	

save "demand_appendix_different_FEs_20210418.dta", replace


*******************
* output
*******************
//Appendix Table 11: demand, different fixed effects 
esttab m1 m2 m3 m4 m5, b(%7.3f) se(%7.3f) stats(N r2_a) br compress nogap   ///
	order(price logshare_cyjinso logshare_cyoins tax energycost1 loghw logw logsize) star(* 0.1 ** 0.05 *** 0.01)
estimates tab m1 m2 m3 m4 m5, b(%7.3f) se(%7.3f) stats(N r2_a) 

esttab m6 m7 m8 m9 m10, b(%7.3f) se(%7.3f) stats(N r2_a) br compress nogap   ///
	order(price logshare_cyjinso logshare_cyoins tax energycost1 loghw logw logsize) star(* 0.1 ** 0.05 *** 0.01)
estimates tab m6 m7 m8 m9 m10, b(%7.3f) se(%7.3f) stats(N r2_a) 

esttab m11 m12 m13, b(%7.3f) se(%7.3f) stats(N r2_a) br compress nogap   ///
	order(price logshare_cyjinso logshare_cyoins tax energycost1 loghw logw logsize) star(* 0.1 ** 0.05 *** 0.01)
estimates tab m11 m12 m13, b(%7.3f) se(%7.3f) stats(N r2_a) 

//elasticity for different fixed effects 
matrix drop A
sum avgela_price_m1, d 
matrix A=(r(p50)\r(mean)\r(sd)\r(p5)\r(p95))
foreach x in m2 m3 m4 m5 m6 m7 m8 m9 m10 m11 m12 m13 {
    sum avgela_price_`x', d
	matrix A=(A,(r(p50)\r(mean)\r(sd)\r(p5)\r(p95)))
}
matrix list A



********************************************************************************
**# Appendix Figure 4 and 5: trend of quality
********************************************************************************
/* Note: in the demand estimation above, the quality does not include country by year fixed
effects, however, when we draw the trend of quality, we need to add the country by year fixed effects back
*/
use "data.dta", clear
* 1.1 MNL_OLS
reghdfe logshare_jcy price tax energycost1 loghw logw logsize  ,   ///
	absorb(i.country2#i.year i.groupmodel2#i.bodytype_harmonized4    ///
	 i.segment2#i.year)  vce(cluster group_mdltrim)
est store ols
matrix beta_ols=e(b)'
matrix list beta_ols 

gen newgamma_jcy_ols=gamma_jcy_ols+fe_cy_ols
gen newwtp_jcy_ols=newgamma_jcy_ols/(-1*beta_ols[1,1])

sum gamma_jcy_ols newgamma_jcy_ols wtp_jcy_ols newwtp_jcy_ols

* 1.2 MNL_IV
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize    ///
	(price=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff),   ///
	absorb(i.country2#i.year i.groupmodel2#i.bodytype_harmonized4    ///
	 i.segment2#i.year) first ffirst cluster(group_mdltrim)
est store mnliv
matrix beta_mnliv=e(b)'
matrix list beta_mnliv

gen newgamma_jcy_mnliv2=gamma_jcy_mnliv2+fe_cy_mnliv2
gen newwtp_jcy_mnliv2=newgamma_jcy_mnliv2/beta_mnliv[1,1]*-1

sum gamma_jcy_mnliv2 newgamma_jcy_mnliv2 wtp_jcy_mnliv2 newwtp_jcy_mnliv2

* 1.3 NL_IV
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(i.country2#i.year i.groupmodel2#i.bodytype_harmonized4    ///
	 i.segment2#i.year) first ffirst cluster(group_mdltrim)
est store nl
matrix beta=e(b)'
matrix list beta

gen newgamma_jcy1t=gamma_jcy1t+fe_cy_nlt
gen newwtp_jcy1t=newgamma_jcy1t/beta[1,1]*-1

sum gamma_jcy1t newgamma_jcy1t wtp_jcy1t newwtp_jcy1t

bysort country year: egen newgamma_cy_ols=mean(newgamma_jcy_ols)
bysort country year: egen newwtp_cy_ols=mean(newwtp_jcy_ols)
bysort country year: egen newgamma_cy_mnliv2=mean(newgamma_jcy_mnliv2)
bysort country year: egen newwtp_cy_mnliv2=mean(newwtp_jcy_mnliv2)
bysort country year: egen newgamma_cy1t=mean(newgamma_jcy1t)
bysort country year: egen newwtp_cy1t=mean(newwtp_jcy1t)

save "demand_appendix_newgamma_20201112.dta", replace

* pattern of quality
duplicates drop country year, force
* value
twoway  (scatter newwtp_cy_ols year, ms(circle) msize(vsmall) connect(l) yaxis(1) xlabel(2005(3)2017)  xline(2009 2012, lpattern(dash) lcolor(gray)))    ///
	(scatter newwtp_cy_mnliv2 year, ms(triangle) msize(vsmall) connect(l) yaxis(2) xlabel(2005(3)2017)   xline(2009 2012, lpattern(dash) lcolor(gray)))    ///
	(scatter newwtp_cy1t year, ms(square) msize(vsmall) connect(l) yaxis(3) xlabel(2005(3)2017)  xline(2009 2012, lpattern(dash) lcolor(gray))) , by(country)     ///
	legend(lab(1 "Multinomial Logit (OLS)") lab(2 "Multinomial Logit (IV)") lab(3 "Nested Logit (IV)"))

* index: baseline is 2005 year
sort country year
foreach x in wtp_cy_ols wtp_cy_mnliv2 wtp_cy1t{
    sort country year
	bysort country: gen index_new`x'=new`x'/new`x'[1]
}
sum index_new*
br country year wtp_cy_ols wtp_cy_mnliv2 wtp_cy1t index_*

twoway  (scatter index_newwtp_cy_ols year, ms(circle) msize(vsmall) connect(l)  xlabel(2005(3)2017)  xline(2009 2012, lpattern(dash) lcolor(gray)) ylabel(1(0.02)1.1))    ///
	(scatter index_newwtp_cy_mnliv2 year, ms(triangle) msize(vsmall) connect(l)  xlabel(2005(3)2017)   xline(2009 2012, lpattern(dash) lcolor(gray)) ylabel(1(0.02)1.1))    ///
	(scatter index_newwtp_cy1t year, ms(square) msize(vsmall) connect(l)  xlabel(2005(3)2017)  xline(2009 2012, lpattern(dash) lcolor(gray)) ylabel(1(0.02)1.1)) , by(country)     ///
	legend(lab(1 "Multinomial Logit (OLS)") lab(2 "Multinomial Logit (IV)") lab(3 "Nested Logit (IV)"))

sum wtp_cy* if year==2005

*******************************
* demand: other specifications of demand estimation
*******************************
use "demand_appendix_newgamma_20201112.dta", clear
sort country year vehicle

* 2.1 NL_s: one-level nest: segment
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize    ///
	(price logshare_cyjins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff),   ///
	absorb(i.country2#i.year i.groupmodel2#i.bodytype_harmonized4    ///
	 i.segment2#i.year) first ffirst cluster(group_mdltrim)
est store nls
matrix beta_nls=e(b)'
matrix list beta_nls

gen newgamma3_jcy_nls2=gamma3_jcy_nls2+fe_cy_nls2
gen newwtp3_jcy_nls2=newgamma3_jcy_nls2/beta_nls[1,1]*-1

sum gamma3_jcy_nls2 newgamma3_jcy_nls2 wtp3_jcy_nls2 newwtp3_jcy_nls2

bysort country year: egen newgamma3_cy_nls2=mean(newgamma3_jcy_nls2)
bysort country year: egen newwtp3_cy_nls2=mean(newwtp3_jcy_nls2)

twoway line wtp3_cy_nls2 newwtp3_cy_nls2 year, by(country)

* 2.2 NL_o: one-level nest: origin
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize    ///
	(price logshare_cyjino=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff   ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_nest_o sum_width_nest_o sum_height_nest_o    ///
	sum_enginecyl_nest_o),   ///
	absorb(i.country2#i.year i.groupmodel2#i.bodytype_harmonized4    ///
	 i.segment2#i.year) first ffirst cluster(group_mdltrim)
est store nlo	
matrix beta_nlo=e(b)'
matrix list beta_nlo

gen newgamma3_jcy_nlo2=gamma3_jcy_nlo2+fe_cy_nlo2
gen newwtp3_jcy_nlo2=newgamma3_jcy_nlo2/beta_nlo[1,1]*-1

sum gamma3_jcy_nlo2 newgamma3_jcy_nlo2 wtp3_jcy_nlo2 newwtp3_jcy_nlo2

bysort country year: egen newgamma3_cy_nlo2=mean(newgamma3_jcy_nlo2)
bysort country year: egen newwtp3_cy_nlo2=mean(newwtp3_jcy_nlo2)

twoway line wtp3_cy_nlo2 newwtp3_cy_nlo2 year, by(country)

* 2.3 NL_so:one-level nest: segment-origin
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize    ///
	(price logshare_cyjinso=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_nest_so sum_width_nest_so sum_height_nest_so    ///
	sum_enginecyl_nest_so),   ///
	absorb(i.country2#i.year i.groupmodel2#i.bodytype_harmonized4    ///
	 i.segment2#i.year) first ffirst cluster(group_mdltrim)
est store nlso	
matrix beta_nlso=e(b)'
matrix list beta_nlso

gen newgamma3_jcy_nlso2=gamma3_jcy_nlso2+fe_cy_nlso2
gen newwtp3_jcy_nlso2=newgamma3_jcy_nlso2/beta_nlso[1,1]*-1
	
sum gamma3_jcy_nlso2 newgamma3_jcy_nlso2 wtp3_jcy_nlso2 newwtp3_jcy_nlso2

bysort country year: egen newgamma3_cy_nlso2=mean(newgamma3_jcy_nlso2)
bysort country year: egen newwtp3_cy_nlso2=mean(newwtp3_jcy_nlso2)

twoway line wtp3_cy_nlso2 newwtp3_cy_nlso2 year, by(country)


* 2.4 model1 (main regression):two-level nest: segment-->origin, ivfe, model-bodytype FE, clustered mdltrim
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(i.country2#i.year i.groupmodel2#i.bodytype_harmonized4    ///
	 i.segment2#i.year) first ffirst cluster(group_mdltrim)
est store m1
matrix beta_m1=e(b)'
matrix list beta_m1

gen newgamma3_jcy_m1t=gamma3_jcy_m1t+fe_cy_m1t
gen newwtp3_jcy_m1t=newgamma3_jcy_m1t/beta_m1[1,1]*-1

sum gamma3_jcy_m1t newgamma3_jcy_m1t wtp3_jcy_m1t newwtp3_jcy_m1t

bysort country year: egen newgamma3_cy_m1t=mean(newgamma3_jcy_m1t)
bysort country year: egen newwtp3_cy_m1t=mean(newwtp3_jcy_m1t)

twoway line wtp3_cy_m1t newwtp3_cy_m1t year, by(country)


* 2.5 model2:two-level nest: segment-->origin, ivfe, model-bodytype-trimlevel FE, clustered mdltrim
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize    ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(i.country2#i.year i.groupmodel2#i.bodytype_harmonized4#i.trimlevel2    ///
	 i.segment2#i.year) first ffirst cluster(group_mdltrim)
est store m2
matrix beta_m2=e(b)'
matrix list beta_m2

gen newgamma3_jcy_m2t=gamma3_jcy_m2t+fe_cy_m2t
gen newwtp3_jcy_m2t=newgamma3_jcy_m2t/beta_m2[1,1]*-1

sum gamma3_jcy_m2t newgamma3_jcy_m2t wtp3_jcy_m2t newwtp3_jcy_m2t

bysort country year: egen newgamma3_cy_m2t=mean(newgamma3_jcy_m2t)
bysort country year: egen newwtp3_cy_m2t=mean(newwtp3_jcy_m2t)

twoway line wtp3_cy_m2t newwtp3_cy_m2t year, by(country)


* 2.6 model4:two-level nest: segment-->origin, ivfe, model-bodytype FE, trimlevel FE, cluster(mdltrim)
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize    ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(i.country2#i.year i.trimlevel2 i.bodytype_harmonized4    ///
	i.groupmodel2#i.bodytype_harmonized4 i.segment2#i.year) first ffirst cluster(group_mdltrim)	
est store m4
matrix beta_m4=e(b)'
matrix list beta_m4

gen newgamma3_jcy_m4t=gamma3_jcy_m4t+fe_cy_m4t
gen newwtp3_jcy_m4t=newgamma3_jcy_m4t/beta_m4[1,1]*-1

sum gamma3_jcy_m4t newgamma3_jcy_m4t wtp3_jcy_m4t newwtp3_jcy_m4t

bysort country year: egen newgamma3_cy_m4t=mean(newgamma3_jcy_m4t)
bysort country year: egen newwtp3_cy_m4t=mean(newwtp3_jcy_m4t)

twoway line wtp3_cy_m4t newwtp3_cy_m4t year, by(country)


* 2.7 model 5: nl_iv, nl for two nesting levels, iv, model-bodytype FE, cluster se, model-year FE
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(i.country2#i.year i.groupmodel2#i.bodytype_harmonized4    ///
	 i.groupmodel2#i.year i.segment2#i.year)    ///
	first ffirst cluster(group_mdltrim)
est store m5
matrix beta_m5=e(b)'
matrix list beta_m5

gen newgamma3_jcy_m5t=gamma3_jcy_m5t+fe_cy_m5t
gen newwtp3_jcy_m5t=newgamma3_jcy_m5t/beta_m5[1,1]*-1

sum gamma3_jcy_m5t newgamma3_jcy_m5t wtp3_jcy_m5t newwtp3_jcy_m5t

bysort country year: egen newgamma3_cy_m5t=mean(newgamma3_jcy_m5t)
bysort country year: egen newwtp3_cy_m5t=mean(newwtp3_jcy_m5t)

twoway line wtp3_cy_m5t newwtp3_cy_m5t year, by(country)


* 2.8 model 6: nl_iv, model-bodytype FE, cluster se, congestion effect
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize lognum_so r_jinso    ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(i.country2#i.year i.groupmodel2#i.bodytype_harmonized4    ///
	 i.segment2#i.year) first ffirst cluster(group_mdltrim)
est store m6
matrix beta_m6=e(b)'
matrix list beta_m6

gen newgamma3_jcy_m6t=gamma3_jcy_m6t+fe_cy_m6t
gen newwtp3_jcy_m6t=newgamma3_jcy_m6t/beta_m6[1,1]*-1

sum gamma3_jcy_m6t newgamma3_jcy_m6t wtp3_jcy_m6t newwtp3_jcy_m6t

bysort country year: egen newgamma3_cy_m6t=mean(newgamma3_jcy_m6t)
bysort country year: egen newwtp3_cy_m6t=mean(newwtp3_jcy_m6t)

twoway line wtp3_cy_m6t newwtp3_cy_m6t year, by(country)

* 2.9 model3: similar to model 1 but with robust se
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(i.country2#i.year i.groupmodel2#i.bodytype_harmonized4    ///
	 i.segment2#i.year) first ffirst robust
est store m3
matrix beta_m3=e(b)'
matrix list beta_m3

gen newgamma3_jcy_m3t=gamma3_jcy_m3t+fe_cy_m3t
gen newwtp3_jcy_m3t=newgamma3_jcy_m3t/beta_m3[1,1]*-1

sum gamma3_jcy_m3t newgamma3_jcy_m3t wtp3_jcy_m3t newwtp3_jcy_m3t

bysort country year: egen newgamma3_cy_m3t=mean(newgamma3_jcy_m3t)
bysort country year: egen newwtp3_cy_m3t=mean(newwtp3_jcy_m3t)

twoway line wtp3_cy_m3t newwtp3_cy_m3t year, by(country)

* 2.10 model 8: allow the preference parameters to vary across time
sort country year vehicle
ivreghdfe logshare_jcy tax i.year#c.energycost1 i.year#c.loghw i.year#c.logw i.year#c.logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(i.country2#i.year i.groupmodel2#i.bodytype_harmonized4     ///
	 i.segment2#i.year) first ffirst cluster(group_mdltrim)
est store m8
matrix beta_m8=e(b)'
matrix list beta_m8

gen newgamma3_jcy_m8t=gamma3_jcy_m8t+fe_cy_m8t
gen newwtp3_jcy_m8t=newgamma3_jcy_m8t/beta_m8[1,1]*-1

sum gamma3_jcy_m8t newgamma3_jcy_m8t wtp3_jcy_m8t newwtp3_jcy_m8t

bysort country year: egen newgamma3_cy_m8t=mean(newgamma3_jcy_m8t)
bysort country year: egen newwtp3_cy_m8t=mean(newwtp3_jcy_m8t)

twoway line wtp3_cy_m8t newwtp3_cy_m8t year, by(country)

save "demand_appendix_newgamma_20201112.dta", replace

************************************
* demand: Dummy for luxury cars
************************************
use "demand_luxury_20200830.dta", clear
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize routlier3_price  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(i.country2#i.year i.groupmodel2#i.bodytype_harmonized4 		///
	 i.segment2#i.year) first ffirst cluster(group_mdltrim)
est store nl_d1
matrix beta_d1=e(b)'
matrix list beta_d1   

gen newgamma3_jcy_d1t=gamma3_jcy_d1t+fe_cy_d1t
gen newwtp3_jcy_d1t=newgamma3_jcy_d1t/beta_d1[1,1]*-1

sum gamma3_jcy_d1t newgamma3_jcy_d1t wtp3_jcy_d1t newwtp3_jcy_d1t

bysort country year: egen newgamma3_cy_d1t=mean(newgamma3_jcy_d1t)
bysort country year: egen newwtp3_cy_d1t=mean(newwtp3_jcy_d1t)

* this looks good
save "demand_luxury_newgamma_20201112.dta", replace

use "demand_luxury_newgamma_20201112.dta", clear
duplicates drop country year, force 
save "demand_luxury_newgamma_countryyear_20201112.dta", replace

************************************
* Output for Appendix Figure 4
************************************
use "demand_appendix_newgamma_20201112.dta", clear

* pattern of quality
duplicates drop country year, force
* value
twoway  (scatter newwtp_cy_ols year, ms(circle) msize(vsmall) connect(l) yaxis(1) xlabel(2005(3)2017)  xline(2009 2012, lpattern(dash) lcolor(gray)))    ///
	(scatter newwtp_cy_mnliv2 year, ms(triangle) msize(vsmall) connect(l) yaxis(2) xlabel(2005(3)2017)   xline(2009 2012, lpattern(dash) lcolor(gray)))    ///
	(scatter newwtp_cy1t year, ms(square) msize(vsmall) connect(l) yaxis(3) xlabel(2005(3)2017)  xline(2009 2012, lpattern(dash) lcolor(gray))) , by(country)     ///
	legend(lab(1 "Multinomial Logit (OLS)") lab(2 "Multinomial Logit (IV)") lab(3 "Nested Logit (IV)"))

* index: baseline is 2005 year
sort country year
foreach x in wtp_cy_ols wtp_cy_mnliv2 wtp_cy1t{
    sort country year
	bysort country: gen index_new`x'=new`x'/new`x'[1]
}
sum index_new*
br country year wtp_cy_ols wtp_cy_mnliv2 wtp_cy1t index_*

twoway  (scatter index_newwtp_cy_ols year, ms(circle) msize(vsmall) connect(l)  xlabel(2005(3)2017)  xline(2009 2012, lpattern(dash) lcolor(gray)) ylabel(1(0.02)1.1))    ///
	(scatter index_newwtp_cy_mnliv2 year, ms(triangle) msize(vsmall) connect(l)  xlabel(2005(3)2017)   xline(2009 2012, lpattern(dash) lcolor(gray)) ylabel(1(0.02)1.1))    ///
	(scatter index_newwtp_cy1t year, ms(square) msize(vsmall) connect(l)  xlabel(2005(3)2017)  xline(2009 2012, lpattern(dash) lcolor(gray)) ylabel(1(0.02)1.1)) , by(country)     ///
	legend(lab(1 "Multinomial Logit (OLS)") lab(2 "Multinomial Logit (IV)") lab(3 "Nested Logit (IV)"))

sum newwtp_cy* if year==2005

************************************
* Output for Appendix Figure 5
************************************
* sub-figure: above
use "demand_appendix_newgamma_20201112.dta", clear
duplicates drop country year, force
sum newwtp3_cy_m1t newwtp3_cy_nls2 newwtp3_cy_nlo2 newwtp3_cy_nlso2
sort country year
foreach x in newwtp3_cy_m1t newwtp3_cy_nls2 newwtp3_cy_nlo2 newwtp3_cy_nlso2{
    sort country year
	bysort country: gen index_`x'=`x'/`x'[1]
}
sum index_new*

twoway  (scatter index_newwtp3_cy_m1t year, ms(circle) msize(vsmall) connect(l) yaxis(1))    ///
	(scatter index_newwtp3_cy_nls2 year, ms(triangle) msize(vsmall) connect(l) yaxis(1))    ///
	(scatter index_newwtp3_cy_nlo2 year, ms(square) msize(vsmall) connect(l) yaxis(1))    ///
	(scatter index_newwtp3_cy_nlso2 year, ms(diamond) msize(vsmall) connect(l) yaxis(1))    ///
	,by(country) legend(lab(1 "Nests: Segment-Origin") lab(2 "One Nest: Segment") lab(3 "One Nest: Origin") lab(4 "One Nest: Segment-Origin"))    ///
	xlabel(2005(3)2017)  xline(2009 2012, lpattern(dash) lcolor(gray))

* sub-figure: bottom
use "demand_appendix_newgamma_20201112.dta", clear
sort country year vehicle
sum fe_cy_m8  
matrix list beta_m8

gen gamma3_jcy_m8t_new=dep_m8-(-.00663473)*tax-(-0.08314)*energycost1-0.819429*loghw-    ///
	1.962571*logw-fe_cy_m1t-(-4.94481)
gen wtp3_jcy_m8t_new=gamma3_jcy_m8t_new/0.08703991
gen newgamma3_jcy_m8t_new=gamma3_jcy_m8t_new+fe_cy_m1t
gen newwtp3_jcy_m8t_new=newgamma3_jcy_m8t_new/0.08703991

sum gamma3_jcy_m8t_new newgamma3_jcy_m8t_new  wtp3_jcy_m8t_new newwtp3_jcy_m8t_new

bysort country year: egen wtp3_cy_m8t_new=mean(wtp3_jcy_m8t_new)
bysort year: egen wtp3_y_m8t_new=mean(wtp3_jcy_m8t_new)
bysort country year: egen newgamma3_cy_m8t_new=mean(newgamma3_jcy_m8t_new)
bysort country year: egen newwtp3_cy_m8t_new=mean(newwtp3_jcy_m8t_new)

twoway line wtp3_cy_m8t_new newwtp3_cy_m8t_new year, by(country)

duplicates drop country year, force
merge 1:1 country year using "demand_luxury_newgamma_countryyear_20201112.dta", keepusing(wtp_cy_d1 wtp_cy_d1t wtp3_cy_d1t newwtp3_cy_d1t)
drop _merge

sort country year
sum newwtp3_cy_m1t newwtp3_cy_m2t newwtp3_cy_m4t newwtp3_cy_m5t newwtp3_cy_m6t newwtp3_cy_d1t newwtp3_cy_m8t_new
foreach x in newwtp3_cy_m1t newwtp3_cy_m2t newwtp3_cy_m4t newwtp3_cy_m5t newwtp3_cy_m6t newwtp3_cy_d1t newwtp3_cy_m8t_new{
    sort country year
	bysort country: gen index_`x'=`x'/`x'[1]
}

sum index_new*

twoway  (scatter index_newwtp3_cy_m1t year, ms(circle) msize(vsmall) connect(l) yaxis(1))    ///
	(scatter index_newwtp3_cy_m2t year, ms(triangle) msize(vsmall) connect(l) yaxis(1))    ///
	(scatter index_newwtp3_cy_m4t year, ms(square) msize(vsmall) connect(l) yaxis(1))    ///
	(scatter index_newwtp3_cy_m5t year, ms(diamond) msize(vsmall) connect(l) yaxis(1))    ///
	(scatter index_newwtp3_cy_m6t year, ms(X) msize(small) connect(l) yaxis(1))    ///
	(scatter index_newwtp3_cy_d1t year, ms(plus) msize(small) connect(l) yaxis(1))    ///
	(scatter index_newwtp3_cy_m8t_new year, ms(V) msize(small) connect(l) yaxis(1))    ///
	, by(country) legend(lab(1 "Model 1") lab(2 "Model 2") lab(3 "Model 3")    ///
	lab(4 "Model 4") lab(5 "Model 5") lab(6 "Model 6") lab(7 "Model 7") )    ///
	xlabel(2005(3)2017)  xline(2009 2012, lpattern(dash) lcolor(gray))




**********************************************************************************************************************************
**#                        Part II: Appendix-----DID Regression
**********************************************************************************************************************************
********************************************************************************
**# Appendix Table 14: Tests for Common Pre-Policy Trends
********************************************************************************
use "data.dta", clear
keep if post==1
foreach x in wtp_jcy1 loghw3 logw2 co2emit logprice{
	reghdfe `x' c.stringency_f_m3#c.year stringency_f_m3   ///
		[fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year) vce(cluster group_mdltrim) 
	est store wf0_`x'
	count if e(sample)==1
	matrix n_wf0_`x'=r(N)
	test c.stringency_f_m3#c.year
	matrix test_wf0_`x'=(r(drop),  r(df_r), r(F), r(df), r(p))'
}

//output
esttab wf0_wtp_jcy1 wf0_loghw3 wf0_logw2 wf0_logprice ,     ///
	b(%7.3f) se(%7.3f) stats(N r2_a) br compress  star(* 0.1 ** 0.05 *** 0.01)

matrix A=(test_wf0_wtp_jcy1, test_wf0_loghw3, test_wf0_logw2, test_wj0_logprice)
matrix list A

matrix B=(n_wf0_wtp_jcy1, n_wf0_loghw3, n_wf0_logw2, n_wj0_logprice)
matrix list B

********************************************************************************
**# Appendix Table 15, 18, 19, 20: Main Robustness Results for Quality
********************************************************************************
use "data.dta", clear

ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(i.country2#i.year i.groupmodel2#i.bodytype_harmonized4    ///
	 i.segment2#i.year) first ffirst cluster(group_mdltrim)
est store nl
matrix beta=e(b)'
matrix list beta
gen snl=e(sample)

//normalize the loghw and logw
drop loghw3 logw2
gen loghw3=loghw/abs(beta[1,1])
gen logw2=logw/abs(beta[1,1])
label var loghw3 "loghw (hp/kg) normalized by the absolute value of price coefficient"
label var logw2 "logw (ton) normalized by the absolute value of price coefficient"

* main regressions
foreach x in wtp_jcy1 loghw3 logw2 logprice{
	reghdfe `x' c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3   ///
		[fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year) vce(cluster group_mdltrim) 
	est store wf0_`x'
	count if e(sample)==1
	matrix n_wf0_`x'=r(N)
	test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
	matrix test_wf0_`x'=(r(drop),  r(df_r), r(F), r(df), r(p))'
}


* model 2: control the pre-trend with change in quality, fuelcons, and horsepower
sum delta_*_j
gen delta_wtp_j=delta_gamma_j/abs(beta[1,1])

foreach x in wtp_jcy1 loghw3 logw2 logprice {
	foreach v in wtp hp fuelcons{
		reghdfe `x' c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3    ///
			c.delta_`v'_j#2006.year c.delta_`v'_j#2007.year c.delta_`v'_j#2008.year c.delta_`v'_j#2009.year    ///
			c.delta_`v'_j#2010.year c.delta_`v'_j#2011.year c.delta_`v'_j#2012.year c.delta_`v'_j#2013.year    ///
			c.delta_`v'_j#2014.year c.delta_`v'_j#2015.year c.delta_`v'_j#2016.year c.delta_`v'_j#2017.year    ///
			[fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year) vce(cluster group_mdltrim)
		est store wf_`x'_`v'
		count if e(sample)==1
		matrix n_wf_`x'_`v'=r(N)
		test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
		matrix test_wf_`x'_`v'=(r(drop),  r(df_r), r(F), r(df), r(p))'
	}
}

	
* model 4: median regression
gen str_post2=stringency_f_m3*cond(post==2,1,0)
gen str_post3=stringency_f_m3*cond(post==3,1,0)

foreach x in wtp_jcy1 loghw3 logw2 logprice str_post2 str_post3 {
	reghdfe `x' [fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year, 		///
		resid(newres_`x') ) vce(cluster group_mdltrim) 
}

foreach x in wtp_jcy1 loghw3 logw2 logprice{
	qreg newres_`x' newres_str_post2 newres_str_post3 [fw=reg], vce(robust)
	est store qr50_`x'
	count if e(sample)==1
	matrix n_qr50_`x'=r(N)
	test newres_str_post2 newres_str_post3
	matrix test_qr50_`x'=(r(drop),  r(df_r), r(F), r(df), r(p))'
}

* model 5: use log of vehicle sales as dependent variable
gen logreg_jcy=log(reg)

reghdfe wtp_jcy1 c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3   ///
		[fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year i.segment2#i.country2#i.year) vce(cluster group_mdltrim) 
gen sm7_wtp_jcy1=e(sample)

foreach x in logreg_jcy {
	reghdfe `x' c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3   ///
		[fw=reg] if sm7_wtp_jcy1==1, absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year) vce(cluster group_mdltrim) 
	est store wf0_`x'
	count if e(sample)==1
	matrix n_wf0_`x'=r(N)
	gen swf0_`x'=e(sample)
	test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
	matrix test_wf0_`x'=(r(drop),  r(df_r), r(F), r(df), r(p))'
}


* Model 6: Dummy for luxury cars
use "demand_luxury_20200830.dta", clear
//demand
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize routlier3_price  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(i.country2#i.year i.groupmodel2#i.bodytype_harmonized4 		///
	 i.segment2#i.year) first ffirst cluster(group_mdltrim)
est store nl_d1
matrix beta_d1=e(b)'
matrix list beta_d1   

//normalize the dependent variables
gen loghw_d1=loghw/abs(beta_d1[1,1])
gen logw_d1=logw/abs(beta_d1[1,1])
egen outlierpost3=group(routlier3_price post)

//main regressions
foreach x in wtp_jcy_d1 loghw_d1 logw_d1 logprice{
	reghdfe `x' c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3   ///
		i.outlierpost3		///
		[fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year) vce(cluster group_mdltrim) 
	est store d1_`x'
	count if e(sample)==1
	matrix n_d1_`x'=r(N)
	test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
	matrix test_d1_`x'=(r(drop),  r(df_r), r(F), r(df), r(p))'
}



*********************
* Output
*********************
//Appendix Table 15: Robustness Checks for quality
esttab wf0_wtp_jcy1 wf0_logreg_jcy wf_wtp_jcy1_wtp wf_wtp_jcy1_hp wf_wtp_jcy1_fuelcons qr50_wtp_jcy1 d1_wtp_jcy_d1,     ///
	b(%7.3f) se(%7.3f) stats(N r2_a) nostar br compress order(*.post#* newres_str_post* stringency_* _cons) keep(*.post#* newres_str_post* stringency_* _cons )
matrix A=(test_wf0_wtp_jcy1, test_wf_wtp_jcy1_wtp, test_wf_wtp_jcy1_hp, test_wf_wtp_jcy1_fuelcons, test_qr50_wtp_jcy1, test_d1_wtp_jcy1_d1)
matrix list A
matrix B=(n_wf0_wtp_jcy1, n_wf_wtp_jcy1_wtp, n_wf_wtp_jcy1_hp, n_wf_wtp_jcy1_fuelcons, n_qr50_wtp_jcy1 n_d1_wtp_jcy_d1)
matrix list B


//Appendix Table 18: Robustness Checks for log horsepower to weight
esttab wf0_loghw3 wf_loghw3_wtp wf_loghw3_hp wf_loghw3_fuelcons  qr50_loghw3 d1_loghw_d1,     ///
	b(%7.3f) se(%7.3f) stats(N r2_a) nostar br compress order(*.post#* newres_str_post* stringency_* _cons) keep(*.post#* newres_str_post* stringency_* _cons )
matrix A=(test_wf0_loghw3, test_wf_loghw3_wtp, test_wf_loghw3_hp, test_wf_loghw3_fuelcons, test_qr50_loghw3, test_d1_loghw_d1)
matrix list A
matrix B=(n_wf0_loghw3, n_wf_loghw3_wtp, n_wf_loghw3_hp, n_wf_loghw3_fuelcons, n_qr50_loghw3, n_d1_loghw_d1)
matrix list B
	
//Appendix Table 19: Robustness Checks for log weight
esttab wf0_logw2 wf_logw2_wtp wf_logw2_hp wf_logw2_fuelcons  qr50_logw2 d1_logw_d1 ,     ///
	b(%7.3f) se(%7.3f) stats(N r2_a) nostar br compress order(*.post#* newres_str_post* stringency_* _cons) keep(*.post#* newres_str_post* stringency_* _cons )
matrix A=(test_wf0_logw2, test_wf_logw2_wtp, test_wf_logw2_hp, test_wf_logw2_fuelcons, test_qr50_logw2, test_d1_logw_d1)
matrix list A
matrix B=(n_wf0_logw2, n_wf_logw2_wtp, n_wf_logw2_hp, n_wf_logw2_fuelcons, n_qr50_logw2, n_d1_logw_d1)
matrix list B
	
//Appendix Table 20: Robustness Checks for logprice
esttab wf0_logprice wf_logprice_wtp wf_logprice_hp wf_logprice_fuelcons qr50_logprice d1_logprice,     ///
	b(%7.4f) se(%7.4f) stats(N r2_a) nostar br compress order(*.post#* newres_str_* stringency_* _cons) keep(*.post#* newres_str_* stringency_* _cons )
matrix A=(test_wf0_logprice, test_wf_logprice_wtp, test_wf_logprice_hp, test_wf_logprice_fuelcons, test_qr50_logprice, test_d1_logprice)
matrix list A
matrix B=(n_wf0_logprice, n_wf_logprice_wtp, n_wf_logprice_hp, n_wf_logprice_fuelcons, n_qr50_logprice, n_d1_logprice)
matrix list B





********************************************************************************
**# Appendix Table 16 and 17: Robustness Results for Quality: Alternative Nesting Structures and specifications
********************************************************************************
use "data.dta", clear
des wtp_jcy*
sum price

* alternative nesting structures and specifications
foreach x in m1 nls nlo nlso m2 m4 m5 m6 m3{
	reghdfe wtp_jcy_`x' c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3   ///
		[fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year) vce(cluster group_mdltrim) 
	est store wtp_`x'
	count if e(sample)==1
	matrix n_wtp_`x'=r(N)
	test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
	matrix test_wtp_`x'=(r(drop),  r(df_r), r(F), r(df), r(p))'
}

* time-variant preference parameters
foreach x in m8{
	reghdfe wtp_jcy_`x' c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3   ///
		[fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year) vce(cluster group_mdltrim) 
	est store wtp_`x'
	count if e(sample)==1
	matrix n_wtp_`x'=r(N)
	test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
	matrix test_wtp_`x'=(r(drop),  r(df_r), r(F), r(df), r(p))'
}

* output
//Appendix Table 16
esttab wtp_m1 wtp_nls wtp_nlo wtp_nlso,  b(%7.3f) se(%7.3f) stats(N r2_a) nostar br compress order(*.post#* stringency_* _cons)
matrix A=(test_wtp_m1, test_wtp_nls, test_wtp_nlo, test_wtp_nlso)
matrix list A
matrix B=(n_wtp_m1, n_wtp_nls, n_wtp_nlo, n_wtp_nlso)
matrix list B
matrix drop A B

//Appendix Table 17
esttab wtp_m1 wtp_m2 wtp_m4 wtp_m5 wtp_m6 wtp_m8 wtp_m3 , b(%7.3f) se(%7.3f) stats(N r2_a) nostar br compress order(*.post#* stringency_* _cons)
matrix A=(test_wtp_m1, test_wtp_m2, test_wtp_m4, test_wtp_m5, test_wtp_m6, test_wtp_m8, test_wtp_m3)
matrix list A
matrix B=(n_wtp_m1, n_wtp_m2, n_wtp_m4, n_wtp_m5, n_wtp_m6, n_wtp_m8, n_wtp_m3)
matrix list B

* dummy for luxury car 
use "demand_luxury_20200830.dta", clear
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize routlier3_price  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(i.country2#i.year i.groupmodel2#i.bodytype_harmonized4 		///
	 i.segment2#i.year) first ffirst cluster(group_mdltrim)
est store nl_d1
matrix beta_d1=e(b)'
matrix list beta_d1   

//normalize the dependent variables
gen loghw_d1=loghw/abs(beta_d1[1,1])
gen logw_d1=logw/abs(beta_d1[1,1])
egen outlierpost3=group(routlier3_price post)

//Quality regression
foreach x in wtp_jcy_d1 {
	reghdfe `x' c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3   ///
		i.outlierpost3		///
		[fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year) vce(cluster group_mdltrim) 
	est store d1_`x'
	count if e(sample)==1
	matrix n_d1_`x'=r(N)
	test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
	matrix test_d1_`x'=(r(drop),  r(df_r), r(F), r(df), r(p))'
}

//output for Appendix Table 17: dummy for luxury cars
esttab d1_wtp_jcy_d1, b(%7.3f) se(%7.3f) stats(N r2_a) br compress nogap star(* 0.1 ** 0.05 *** 0.01)
matrix A=(test_d1_wtp_jcy_d1)
matrix list A
matrix B=(n_d1_wtp_jcy_d1)
matrix list B


********************************************************************************
**# Appendix Table 21-24: DID regressions, different fixed effects
********************************************************************************
use "data.dta", clear

ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(i.country2#i.year i.groupmodel2#i.bodytype_harmonized4    ///
	 i.segment2#i.year) first ffirst cluster(group_mdltrim)
est store nl
matrix beta=e(b)'
matrix list beta
gen snl=e(sample)

//normalize the loghw and logw
drop loghw3 logw2
gen loghw3=loghw/abs(beta[1,1])
gen logw2=logw/abs(beta[1,1])
label var loghw3 "loghw (hp/kg) normalized by the absolute value of price coefficient"
label var logw2 "logw (ton) normalized by the absolute value of price coefficient"

* main regressions
foreach x in wtp_jcy1 loghw3 logw2 logprice{
	reghdfe `x' c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3   ///
		[fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year) vce(cluster group_mdltrim) 
	est store wf0_`x'
	count if e(sample)==1
	matrix n_wf0_`x'=r(N)
	test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
	matrix test_wf0_`x'=(r(drop),  r(df_r), r(F), r(df), r(p))'
}

//model-trim FE
foreach x in wtp_jcy1 loghw3 logw2 logprice{
	reghdfe `x' c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3   ///
		[fw=reg], absorb(group_mdltrim bodytype_harmonized2 numberofdoors     ///
			drivenwheels2 transtype fuelcat enginecyl numberofgears    ///
       i.country2 i.year i.country2#i.year i.segment2#i.year ) vce(cluster group_mdltrim)
	est store m2_`x'
	count if e(sample)==1
	matrix n_m2_`x'=r(N)
	test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
	matrix test_m2_`x'=(r(drop),  r(df_r), r(F), r(df), r(p))'
}


//model-trim-bodytype
egen mt_b=group(group_mdltrim bodytype_harmonized)
foreach x in wtp_jcy1 loghw3 logw2 logprice{
	reghdfe `x' c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3   ///
			[fw=reg], absorb(mt_b numberofdoors     ///
			drivenwheels2 transtype fuelcat enginecyl numberofgears    ///
			i.country2 i.year i.country2#i.year i.segment2#i.year ) vce(cluster group_mdltrim)
	est store m2_`x'_b2
	count if e(sample)==1
	matrix n_m2_`x'_b2=r(N)
	test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
	matrix test_m2_`x'_b2=(r(drop),  r(df_r), r(F), r(df), r(p))'
}


//model-trim-bodytype-fueltype
egen mt_bft=group(group_mdltrim bodytype_harmonized fuelcat)
foreach x in wtp_jcy1 loghw3 logw2 logprice{
	reghdfe `x' c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3   ///
			[fw=reg], absorb(mt_bft numberofdoors     ///
			drivenwheels2 transtype enginecyl numberofgears    ///
			i.country2 i.year i.country2#i.year i.segment2#i.year ) vce(cluster group_mdltrim)
	est store m2_`x'_bft2
	count if e(sample)==1
	matrix n_m2_`x'_bft2=r(N)
	test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
	matrix test_m2_`x'_bft2=(r(drop),  r(df_r), r(F), r(df), r(p))'
}


//model-trim-bodytype-fueltype-segment
egen mt_bfts=group(group_mdltrim bodytype_harmonized fuelcat segment)
foreach x in wtp_jcy1 loghw3 logw2 logprice {
	reghdfe `x' c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3   ///
			[fw=reg], absorb(mt_bfts numberofdoors     ///
			drivenwheels2 transtype enginecyl numberofgears    ///
			i.country2 i.year i.country2#i.year i.segment2#i.year ) vce(cluster group_mdltrim)
	est store m2_`x'_bfts2
	count if e(sample)==1
	matrix n_m2_`x'_bfts2=r(N)
	test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
	matrix test_m2_`x'_bfts2=(r(drop),  r(df_r), r(F), r(df), r(p))'
}


//model-trim-bodytype-fueltype-segment-transtype
egen mt_bftst=group(group_mdltrim bodytype_harmonized fuelcat segment transtype)
foreach x in wtp_jcy1 loghw3 logw2 logprice {
	reghdfe `x' c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3   ///
			[fw=reg], absorb(mt_bftst numberofdoors     ///
			drivenwheels2 enginecyl numberofgears    ///
			i.country2 i.year i.country2#i.year i.segment2#i.year ) vce(cluster group_mdltrim)
	est store m2_`x'_bftst2
	count if e(sample)==1
	matrix n_m2_`x'_bftst2=r(N)
	test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
	matrix test_m2_`x'_bftst2=(r(drop),  r(df_r), r(F), r(df), r(p))'
}

//model-trim-bodytype-fueltype-segment-transtype-drivenwheels
egen mt_bftstd=group(group_mdltrim bodytype_harmonized fuelcat segment transtype drivenwheels)
foreach x in wtp_jcy1 loghw3 logw2 logprice {
	reghdfe `x' c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3   ///
			[fw=reg], absorb(mt_bftstd numberofdoors     ///
			enginecyl numberofgears    ///
			i.country2 i.year i.country2#i.year i.segment2#i.year ) vce(cluster group_mdltrim)
	est store m2_`x'_bftstd2
	count if e(sample)==1
	matrix n_m2_`x'_bftstd2=r(N)
	test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
	matrix test_m2_`x'_bftstd2=(r(drop),  r(df_r), r(F), r(df), r(p))'
}

//model-trim-bodytype-fueltype-segment-transtype-drivenwheels-numberofdoors
egen mt_bftstdn=group(group_mdltrim bodytype_harmonized fuelcat segment transtype drivenwheels numberofdoors)
foreach x in wtp_jcy1 loghw3 logw2 logprice {
	reghdfe `x' c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3   ///
			[fw=reg], absorb(mt_bftstdn      ///
			enginecyl numberofgears    ///
			i.country2 i.year i.country2#i.year i.segment2#i.year ) vce(cluster group_mdltrim)
	est store m2_`x'_bftstdn2
	count if e(sample)==1
	matrix n_m2_`x'_bftstdn2=r(N)
	test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
	matrix test_m2_`x'_bftstdn2=(r(drop),  r(df_r), r(F), r(df), r(p))'
}

//model-trim-bodytype-fueltype-segment-transtype-drivenwheels-numberofdoors-enginecyl
egen mt_bftstdne=group(group_mdltrim bodytype_harmonized fuelcat segment transtype drivenwheels numberofdoors enginecyl)
foreach x in wtp_jcy1 loghw3 logw2 logprice{
	reghdfe `x' c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3   ///
			[fw=reg], absorb(mt_bftstdne numberofgears i.country2 i.year i.country2#i.year i.segment2#i.year ) vce(cluster group_mdltrim)
	est store m2_`x'_bftstdne2
	count if e(sample)==1
	matrix n_m2_`x'_bftstdne2=r(N)
	test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
	matrix test_m2_`x'_bftstdne2=(r(drop),  r(df_r), r(F), r(df), r(p))'
}

**********************
* Output:
**********************
//Appendix Table 21: Quality
esttab m2_wtp_jcy1 m2_wtp_jcy1_b2 m2_wtp_jcy1_bft2 m2_wtp_jcy1_bfts2 m2_wtp_jcy1_bftst2 m2_wtp_jcy1_bftstd2 m2_wtp_jcy1_bftstdn2 m2_wtp_jcy1_bftstdne2 wf0_wtp_jcy1,  b(%7.3f) se(%7.3f) stats(N r2_a) br compress order(*.post#* stringency_* _cons) nostar
estimates tab m2_wtp_jcy1 m2_wtp_jcy1_b2 m2_wtp_jcy1_bft2 m2_wtp_jcy1_bfts2 m2_wtp_jcy1_bftst2 m2_wtp_jcy1_bftstd2 m2_wtp_jcy1_bftstdn2 m2_wtp_jcy1_bftstdne2 wf0_wtp_jcy1,  b(%7.3f) se(%7.3f) stats(N r2_a)
matrix A=(test_m2_wtp_jcy1, test_m2_wtp_jcy1_b2, test_m2_wtp_jcy1_bft2, test_m2_wtp_jcy1_bfts2, 		///
	test_m2_wtp_jcy1_bftst2, test_m2_wtp_jcy1_bftstd2, test_m2_wtp_jcy1_bftstdn2, test_m2_wtp_jcy1_bftstdne2, test_wf0_wtp_jcy1 )
matrix list A
matrix B=(n_m2_wtp_jcy1, n_m2_wtp_jcy1_b2, n_m2_wtp_jcy1_bft2, n_m2_wtp_jcy1_bfts2, 		///
	n_m2_wtp_jcy1_bftst2, n_m2_wtp_jcy1_bftstd2, n_m2_wtp_jcy1_bftstdn2, n_m2_wtp_jcy1_bftstdne2, n_wf0_wtp_jcy1 )
matrix list B
	
//Appendix Table 22: log horsepower / weight
esttab m2_loghw3 m2_loghw3_b2 m2_loghw3_bft2 m2_loghw3_bfts2 m2_loghw3_bftst2 m2_loghw3_bftstd2 m2_loghw3_bftstdn2 m2_loghw3_bftstdne2 wf0_loghw3,  b(%7.3f) se(%7.3f) stats(N r2_a) br compress order(*.post#* stringency_* _cons) nostar
estimates tab m2_loghw3 m2_loghw3_b2 m2_loghw3_bft2 m2_loghw3_bfts2 m2_loghw3_bftst2 m2_loghw3_bftstd2 m2_loghw3_bftstdn2 m2_loghw3_bftstdne2 wf0_loghw3,  b(%7.3f) se(%7.3f) stats(N r2_a) 
matrix A=(test_m2_loghw3, test_m2_loghw3_b2, test_m2_loghw3_bft2, test_m2_loghw3_bfts2, 		///
	test_m2_loghw3_bftst2, test_m2_loghw3_bftstd2, test_m2_loghw3_bftstdn2, test_m2_loghw3_bftstdne2, test_wf0_loghw3 )
matrix list A
matrix B=(n_m2_loghw3, n_m2_loghw3_b2, n_m2_loghw3_bft2, n_m2_loghw3_bfts2, 		///
	n_m2_loghw3_bftst2, n_m2_loghw3_bftstd2, n_m2_loghw3_bftstdn2, n_m2_loghw3_bftstdne2, n_wf0_loghw3 )
matrix list B

//Appendix Table 23: log weight
esttab m2_logw2 m2_logw2_b2 m2_logw2_bft2 m2_logw2_bfts2 m2_logw2_bftst2 m2_logw2_bftstd2 m2_logw2_bftstdn2 m2_logw2_bftstdne2 wf0_logw2,  b(%7.3f) se(%7.3f) stats(N r2_a) br compress order(*.post#* stringency_* _cons) nostar
estimates tab m2_logw2 m2_logw2_b2 m2_logw2_bft2 m2_logw2_bfts2 m2_logw2_bftst2 m2_logw2_bftstd2 m2_logw2_bftstdn2 m2_logw2_bftstdne2 wf0_logw2,  b(%7.3f) se(%7.3f) stats(N r2_a) 
matrix A=(test_m2_logw2, test_m2_logw2_b2, test_m2_logw2_bft2, test_m2_logw2_bfts2, 		///
	test_m2_logw2_bftst2, test_m2_logw2_bftstd2, test_m2_logw2_bftstdn2, test_m2_logw2_bftstdne2, test_wf0_logw2 )
matrix list A
matrix B=(n_m2_logw2, n_m2_logw2_b2, n_m2_logw2_bft2, n_m2_logw2_bfts2, 		///
	n_m2_logw2_bftst2, n_m2_logw2_bftstd2, n_m2_logw2_bftstdn2, n_m2_logw2_bftstdne2, n_wf0_logw2 )
matrix list B

//Appendix Table 24: log price
esttab m2_logprice m2_logprice_b2 m2_logprice_bft2 m2_logprice_bfts2 m2_logprice_bftst2 m2_logprice_bftstd2 m2_logprice_bftstdn2 m2_logprice_bftstdne2 wf0_logprice,  b(%7.3f) se(%7.3f) stats(N r2_a) br compress order(*.post#* stringency_* _cons) nostar
estimates tab m2_logprice m2_logprice_b2 m2_logprice_bft2 m2_logprice_bfts2 m2_logprice_bftst2 m2_logprice_bftstd2 m2_logprice_bftstdn2 m2_logprice_bftstdne2 wf0_logprice,  b(%7.3f) se(%7.3f) stats(N r2_a) 
matrix A=(test_m2_logprice, test_m2_logprice_b2, test_m2_logprice_bft2, test_m2_logprice_bfts2, 		///
	test_m2_logprice_bftst2, test_m2_logprice_bftstd2, test_m2_logprice_bftstdn2, test_m2_logprice_bftstdne2, test_wf0_logprice )
matrix list A
matrix B=(n_m2_logprice, n_m2_logprice_b2, n_m2_logprice_bft2, n_m2_logprice_bfts2, 		///
	n_m2_logprice_bftst2, n_m2_logprice_bftstd2, n_m2_logprice_bftstdn2, n_m2_logprice_bftstdne2, n_wf0_logprice )
matrix list B


********************************************************************************
**# Appendix Figure 6: percentage Change in a Firm's 2008-2011 Registrations
********************************************************************************
use "data.dta", clear

ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(i.country2#i.year i.groupmodel2#i.bodytype_harmonized4    ///
	 i.segment2#i.year) first ffirst cluster(group_mdltrim)
est store nl
matrix beta=e(b)'
matrix list beta
gen snl=e(sample)

//normalize the loghw and logw
drop loghw3 logw2
gen loghw3=loghw/abs(beta[1,1])
gen logw2=logw/abs(beta[1,1])
label var loghw3 "loghw (hp/kg) normalized by the absolute value of price coefficient"
label var logw2 "logw (ton) normalized by the absolute value of price coefficient"

unique stringency_f_m3, by(company_harmonized3)
//firm's change in total sales from 2007-2011 
bysort company_harmonized3 year: egen reg_fy=total(reg), missing
duplicates drop company_harmonized3 year, force	
sort company_harmonized3 year
by company_harmonized3: gen changereg_f0711=reg_fy[7]-reg_fy[3]
br company_harmonized3 year reg_fy changereg_f0711 stringency_f_m3
sort company_harmonized3 year
by company_harmonized3: gen perchangereg_f0711=(reg_fy[7]-reg_fy[3])/reg_fy[3]

bysort year: egen reg_y=total(reg_fy), missing
bysort company_harmonized3 year: gen share_fy=reg_fy/reg_y
br company_harmonized3 year reg_fy reg_y share_fy
sort company_harmonized3 year
by company_harmonized3: gen changeshare_f0711=share_fy[7]-share_fy[3]
sort company_harmonized3 year
by company_harmonized3: gen perchangeshare_f0711=(share_fy[7]-share_fy[3])/share_fy[3]
br company_harmonized3 year share_fy changeshare_f0711 perchangeshare_f0711 stringency_f_m3

//scatter plot
duplicates drop company_harmonized3, force
twoway (scatter changereg_f0711 stringency_f_m3, mlabposition(2) mlabel(company_harmonized3) msize(small)) (lfit changereg_f0711 stringency_f_m3, range(0.1 0.3)), xtitle("firm-level stringency") ytitle("2007-2011 firm's change in total sales")

twoway (scatter perchangereg_f0711 stringency_f_m3, mlabposition(3) mlabel(company_harmonized3) msize(small)) (lfit perchangereg_f0711 stringency_f_m3, range(0.1 0.3)), xtitle("firm-level stringency") ytitle("2007-2011 firm's percentage change in total sales") legend(off)


********************************************************************************
**# Appendix Table 25: Quality regression,Different Divisions of Time Periods
********************************************************************************
use "data.dta", clear
//demand
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(i.country2#i.year i.groupmodel2#i.bodytype_harmonized4    ///
	 i.segment2#i.year) first ffirst cluster(group_mdltrim)
est store nl
matrix beta=e(b)'
matrix list beta
gen snl=e(sample)

//normalize the loghw and logw
drop loghw3 logw2
gen loghw3=loghw/abs(beta[1,1])
gen logw2=logw/abs(beta[1,1])
label var loghw3 "loghw (hp/kg) normalized by the absolute value of price coefficient"
label var logw2 "logw (ton) normalized by the absolute value of price coefficient"
save "demand_appendix_20210114.dta", replace

* (1) main regressions using firm level stringnecy
foreach x in wtp_jcy1 loghw3 logw2 logprice {
	reghdfe `x' c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3   ///
		[fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year) vce(cluster group_mdltrim) 
	est store wf0_`x'
	count if e(sample)==1
	matrix n_wf0_`x'=r(N)
	test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
	matrix test_wf0_`x'=(r(drop),  r(df_r), r(F), r(df), r(p))'
}

* output
esttab wf0_wtp_jcy1 wf0_loghw3 wf0_logw2 wf0_logprice,     ///
	b(%7.3f) se(%7.3f) stats(N r2_a) br compress order(*.post#* stringency_* _cons) star(* 0.1 ** 0.05 *** 0.01) nogap keep(*.post#*)
estimates tab wf0_wtp_jcy1 wf0_loghw3 wf0_logw2, b(%7.3f) se(%7.3f) stats(N r2_a)

matrix A=(test_wf0_wtp_jcy1, test_wf0_loghw3, test_wf0_logw2, test_wf0_logprice)
matrix list A

matrix B=(n_wf0_wtp_jcy1, n_wf0_loghw3, n_wf0_logw2, n_wf0_logprice)
matrix list B


* (2) 2005-2008, 2009-2017
gen post2008=(year>=2009)
* main regressions using firm level stringnecy
foreach x in wtp_jcy1 loghw3 logw2 logprice {
	reghdfe `x' c.stringency_f_m3#post2008 stringency_f_m3 ///
		[fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year) vce(cluster group_mdltrim) 
	est store wf0_`x'
	count if e(sample)==1
	matrix n_wf0_`x'=r(N)
}

* output
esttab wf0_wtp_jcy1 wf0_loghw3 wf0_logw2 wf0_logprice,     ///
	b(%7.3f) se(%7.3f) stats(N r2_a) br compress order(*.post2008#* stringency_* _cons) star(* 0.1 ** 0.05 *** 0.01) nogap keep(1.post2008#* _cons)

matrix B=(n_wf0_wtp_jcy1, n_wf0_loghw3, n_wf0_logw2, n_wf0_logprice)
matrix list B

reghdfe wtp_jcy1 c.stringency_f_m3#post2008 stringency_f_m3 ///
		[fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year) vce(cluster group_mdltrim) 
test c.stringency_f_m3#1.post2008

matrix test_wf0_wtp_jcy1=(r(drop),  r(df_r), r(F), r(df), r(p))'
matrix list test_wf0_wtp_jcy1


* (3) 2005-2007, 2008-2017
use "demand_appendix_20210114.dta", clear
des post*
gen post2007=(year>=2008)

* main regressions using firm level stringnecy
foreach x in wtp_jcy1 loghw3 logw2 logprice {
	reghdfe `x' c.stringency_f_m3#post2007 stringency_f_m3 ///
		[fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year) vce(cluster group_mdltrim) 
	est store wf0_`x'
	count if e(sample)==1
	matrix n_wf0_`x'=r(N)
}

* output
esttab wf0_wtp_jcy1 wf0_loghw3 wf0_logw2 wf0_logprice,     ///
	b(%7.3f) se(%7.3f) stats(N r2_a) br compress order(*.post2007#* stringency_* _cons) star(* 0.1 ** 0.05 *** 0.01) nogap keep(1.post2007#* _cons)

matrix B=(n_wf0_wtp_jcy1, n_wf0_loghw3, n_wf0_logw2, n_wf0_logprice, n_wj0_logprice)
matrix list B

reghdfe wtp_jcy1 c.stringency_f_m3#post2007 stringency_f_m3 ///
		[fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year) vce(cluster group_mdltrim) 
test c.stringency_f_m3#1.post2007
matrix test_wf0_wtp_jcy1=(r(drop),  r(df_r), r(F), r(df), r(p))'
matrix list test_wf0_wtp_jcy1

* (4) 2005-2007, 2008-2011, 2012-2017
use "demand_appendix_20210114.dta", clear
des post*
gen post0712=1 if year<=2007
replace post0712=2 if year>=2008 & year<=2011
replace post0712=3 if year>=2012
tab year post0712

* Main regression: firm-level stringency, drop 1% uppper side outliers of price, segment-year fixed effects
* main regressions using firm level stringnecy
foreach x in wtp_jcy1 loghw3 logw2 logprice {
	reghdfe `x' c.stringency_f_m3#2.post0712 c.stringency_f_m3#3.post0712 stringency_f_m3   ///
		[fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year) vce(cluster group_mdltrim) 
	est store wf0_`x'
	count if e(sample)==1
	matrix n_wf0_`x'=r(N)
	test c.stringency_f_m3#2.post0712 c.stringency_f_m3#3.post0712
	matrix test_wf0_`x'=(r(drop),  r(df_r), r(F), r(df), r(p))'
}

* output
esttab wf0_wtp_jcy1 wf0_loghw3 wf0_logw2 wf0_logprice,     ///
	b(%7.3f) se(%7.3f) stats(N r2_a) br compress order(*.post0712#* stringency_* _cons) star(* 0.1 ** 0.05 *** 0.01) nogap keep(*.post0712#*)

matrix A=(test_wf0_wtp_jcy1, test_wf0_loghw3, test_wf0_logw2, test_wf0_logprice)
matrix list A

matrix B=(n_wf0_wtp_jcy1, n_wf0_loghw3, n_wf0_logw2, n_wf0_logprice)
matrix list B


* (5)) 2005-2009, 2010-2017
use "demand_appendix_20210114.dta", clear
des post*
gen post2010=(year>=2010)

* main regressions using firm level stringnecy
foreach x in wtp_jcy1 loghw3 logw2 logprice {
	reghdfe `x' c.stringency_f_m3#post2010 stringency_f_m3 ///
		[fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year) vce(cluster group_mdltrim) 
	est store wf0_`x'
	count if e(sample)==1
	matrix n_wf0_`x'=r(N)
}

* output
esttab wf0_wtp_jcy1 wf0_loghw3 wf0_logw2 wf0_logprice,     ///
	b(%7.3f) se(%7.3f) stats(N r2_a) br compress order(*.post2010#* stringency_* _cons) star(* 0.1 ** 0.05 *** 0.01) nogap keep(1.post2010#* _cons)

matrix B=(n_wf0_wtp_jcy1, n_wf0_loghw3, n_wf0_logw2, n_wf0_logprice)
matrix list B

reghdfe wtp_jcy1 c.stringency_f_m3#post2010 stringency_f_m3 ///
		[fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year) vce(cluster group_mdltrim) 
test c.stringency_f_m3#1.post2010 
matrix test_wf0_wtp_jcy1=(r(drop),  r(df_r), r(F), r(df), r(p))'
matrix list test_wf0_wtp_jcy1


* (6))2005-2009, 2010-2011, 2012-2017
use "demand_appendix_20210114.dta", clear
des post*
gen post1012=1 if year<=2009
replace post1012=2 if year>=2010 & year<=2011
replace post1012=3 if year>=2012
tab year post1012

* Main regression: firm-level stringency, drop 1% uppper side outliers of price, segment-year fixed effects
* main regressions using firm level stringnecy
foreach x in wtp_jcy1 loghw3 logw2 logprice {
	reghdfe `x' c.stringency_f_m3#2.post1012 c.stringency_f_m3#3.post1012 stringency_f_m3   ///
		[fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year) vce(cluster group_mdltrim) 
	est store wf0_`x'
	count if e(sample)==1
	matrix n_wf0_`x'=r(N)
	test c.stringency_f_m3#2.post1012 c.stringency_f_m3#3.post1012
	matrix test_wf0_`x'=(r(drop),  r(df_r), r(F), r(df), r(p))'
}

* output
esttab wf0_wtp_jcy1 wf0_loghw3 wf0_logw2 wf0_logprice,     ///
	b(%7.3f) se(%7.3f) stats(N r2_a) br compress order(*.post1012#* stringency_* _cons) star(* 0.1 ** 0.05 *** 0.01) nogap keep(*.post1012#*)

matrix A=(test_wf0_wtp_jcy1, test_wf0_loghw3, test_wf0_logw2, test_wf0_logprice)
matrix list A

matrix B=(n_wf0_wtp_jcy1, n_wf0_loghw3, n_wf0_logw2, n_wf0_logprice)
matrix list B

* (7) 2005-2008, 2009-2012, 2013-2017
use "demand_appendix_20210114.dta", clear
des post*

gen post0913=1 if year<=2008
replace post0913=2 if year>=2009 & year<=2012
replace post0913=3 if year>=2013 & year<=2017
tab year post0913

* Main regression: firm-level stringency, drop 1% uppper side outliers of price, segment-year fixed effects
* main regressions using firm level stringnecy
foreach x in wtp_jcy1 loghw3 logw2 logprice {
	reghdfe `x' c.stringency_f_m3#2.post0913 c.stringency_f_m3#3.post0913 stringency_f_m3   ///
		[fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year) vce(cluster group_mdltrim) 
	est store wf0_`x'
	count if e(sample)==1
	matrix n_wf0_`x'=r(N)
	test c.stringency_f_m3#2.post0913 c.stringency_f_m3#3.post0913
	matrix test_wf0_`x'=(r(drop),  r(df_r), r(F), r(df), r(p))'
}

* output
esttab wf0_wtp_jcy1 wf0_loghw3 wf0_logw2 wf0_logprice,     ///
	b(%7.3f) se(%7.3f) stats(N r2_a) br compress order(*.post0913#* stringency_* _cons) star(* 0.1 ** 0.05 *** 0.01) nogap keep(*.post0913#*)

matrix A=(test_wf0_wtp_jcy1, test_wf0_loghw3, test_wf0_logw2, test_wf0_logprice)
matrix list A

matrix B=(n_wf0_wtp_jcy1, n_wf0_loghw3, n_wf0_logw2, n_wf0_logprice)
matrix list B

* (8) 2005-2008, 2009-2013, 2014-2017
use "demand_appendix_20210114.dta", clear
des post*
gen post0914=1 if year<=2008
replace post0914=2 if year>=2009 & year<=2013
replace post0914=3 if year>=2014 & year<=2017
tab year post0914

* Main regression: firm-level stringency, drop 1% uppper side outliers of price, segment-year fixed effects
* main regressions using firm level stringnecy
foreach x in wtp_jcy1 loghw3 logw2 logprice {
	reghdfe `x' c.stringency_f_m3#2.post0914 c.stringency_f_m3#3.post0914 stringency_f_m3   ///
		[fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year) vce(cluster group_mdltrim) 
	est store wf0_`x'
	count if e(sample)==1
	matrix n_wf0_`x'=r(N)
	test c.stringency_f_m3#2.post0914 c.stringency_f_m3#3.post0914
	matrix test_wf0_`x'=(r(drop),  r(df_r), r(F), r(df), r(p))'
}

* output
esttab wf0_wtp_jcy1 wf0_loghw3 wf0_logw2 wf0_logprice,     ///
	b(%7.3f) se(%7.3f) stats(N r2_a) br compress order(*.post0914#* stringency_* _cons) star(* 0.1 ** 0.05 *** 0.01) nogap keep(*.post0914#*)

matrix A=(test_wf0_wtp_jcy1, test_wf0_loghw3, test_wf0_logw2, test_wf0_logprice)
matrix list A

matrix B=(n_wf0_wtp_jcy1, n_wf0_loghw3, n_wf0_logw2, n_wf0_logprice)
matrix list B



********************************************************************************
**# Appendix Table 26: Quality regression, Controlling for Scrappage Schemes
********************************************************************************
* 1. define variables for scrapping scheme
use "demand_appendix_20210114.dta", clear
//near the qualification for those targeted scrapping scheme
foreach x in 5 10 15 20{
	gen scrap`x'=1 if country=="France" & year==2008 & co2emit<=160 & co2emit>=160*(1-`x'*0.01)
	replace scrap`x'=1 if country=="France" & year==2009 & co2emit<=160 & co2emit>=160*(1-`x'*0.01)
	replace scrap`x'=1 if country=="France" & year==2010 & co2emit<=155 & co2emit>=155*(1-`x'*0.01)
	replace scrap`x'=1 if country=="France" & year==2011 & co2emit<=150 & co2emit>=150*(1-`x'*0.01)
	replace scrap`x'=1 if country=="Italy" & year==2007 & ((co2emit<=140 & co2emit>=140*(1-`x'*0.01) & fuelcat_harmonized2==3) | (co2emit<=130 & co2emit>=130*(1-`x'*0.01) & fuelcat_harmonized2==1))
	replace scrap`x'=1 if country=="Italy" & year==2008 & co2emit<=130 & co2emit>=130*(1-`x'*0.01)
	replace scrap`x'=1 if country=="Italy" & year==2009 & ((co2emit<=140 & co2emit>=140*(1-`x'*0.01) & fuelcat_harmonized2==3) | (co2emit<=130 & co2emit>=130*(1-`x'*0.01) & fuelcat_harmonized2==1) | fuelcat_harmonized2==2)
	replace scrap`x'=1 if country=="Spain" & year==2007 & engineccm<=2500 & engineccm>=2500*(1-`x'*0.01)
	replace scrap`x'=1 if country=="Spain" & year>=2008 & year<=2010 & co2emit<=120 & co2emit>=120*(1-`x'*0.01) & price<=30
	replace scrap`x'=0 if scrap`x'==.
}

* 2. regressions with dummies for scrapping scheme 
* main regressions using firm level stringnecy
foreach y in 5 10 15 20{
	foreach x in wtp_jcy1{
		reghdfe `x' c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3 scrap`y'    ///
			[fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year) vce(cluster group_mdltrim) 
		est store wf0_`x'_`y'
		count if e(sample)==1
		matrix n_wf0_`x'_`y'=r(N)
		test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
		matrix test_wf0_`x'_`y'=(r(drop),  r(df_r), r(F), r(df), r(p))'
	}
}

* 3. dummy for all eligible vehicles for targeted scrapping scheme
gen scrap=1 if country=="France" & year==2008 & co2emit<=160 
replace scrap=1 if country=="France" & year==2009 & co2emit<=160 
replace scrap=1 if country=="France" & year==2010 & co2emit<=155 
replace scrap=1 if country=="France" & year==2011 & co2emit<=150 
replace scrap=1 if country=="Italy" & year==2007 & ((co2emit<=140 & fuelcat_harmonized2==3) | (co2emit<=130 & fuelcat_harmonized2==1))
replace scrap=1 if country=="Italy" & year==2008 & co2emit<=130 
replace scrap=1 if country=="Italy" & year==2009 & ((co2emit<=140 & fuelcat_harmonized2==3) | (co2emit<=130 & fuelcat_harmonized2==1) | fuelcat_harmonized2==2)
replace scrap=1 if country=="Spain" & year==2007 & engineccm<=2500
replace scrap=1 if country=="Spain" & year>=2008 & year<=2010 & co2emit<=120 & price<=30
replace scrap=0 if scrap==.

* regressions with dummies for scrapping scheme 
foreach x in wtp_jcy1 loghw3 logw2 logprice {
	reghdfe `x' c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3 scrap    ///
			[fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year) vce(cluster group_mdltrim) 
	est store wf0_`x'
	count if e(sample)==1
	matrix n_wf0_`x'=r(N)
	test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
	matrix test_wf0_`x'=(r(drop),  r(df_r), r(F), r(df), r(p))'
}
* output
esttab wf0_wtp_jcy1 wf0_wtp_jcy1_5 wf0_wtp_jcy1_10 wf0_wtp_jcy1_15 wf0_wtp_jcy1_20,     ///
	b(%7.3f) se(%7.3f) stats(N r2_a) br compress order(*.post#* stringency_* _cons) star(* 0.1 ** 0.05 *** 0.01) nogap keep(*.post#* scrap5)


matrix A=(test_wf0_wtp_jcy1, test_wf0_wtp_jcy1_5, test_wf0_wtp_jcy1_10, test_wf0_wtp_jcy1_15, test_wf0_wtp_jcy1_20)
matrix list A

matrix B=(n_wf0_wtp_jcy1, n_wf0_wtp_jcy1_5, n_wf0_wtp_jcy1_10, n_wf0_wtp_jcy1_15, n_wf0_wtp_jcy1_20)
matrix list B



********************************************************************************
**# Appendix Table 27: Quality regression, Different fixed effects for demand
********************************************************************************
* 1. different fixed effects for the demand estimation
use "data.dta", clear

* 1.1 model1 (main regression):two-level nest: segment-->origin, ivfe, model-bodytype FE, clustered mdltrim
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m1=i.country2#i.year fe_mb_m1=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_m1=i.segment2#i.year, resid(resid_m1)) first ffirst cluster(group_mdltrim)
est store m1
matrix beta_m1=e(b)'
matrix list beta_m1

gen gamma_jcy_m1=logshare_jcy-beta_m1[1,1]*price-beta_m1[2,1]*logshare_cyjinso-beta_m1[3,1]*logshare_cyoins    ///
	-beta_m1[4,1]*tax-beta_m1[5,1]*energycost1-beta_m1[6,1]*loghw-beta_m1[7,1]*logw-fe_cy_m1
gen wtp_jcy_m1=gamma_jcy_m1/beta_m1[1,1]*-1

//elasticity and wtp for attributes
gen elastic_price_m1=((1-share_cyjinso)/(1-beta_m1[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m1[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m1[1,1]*price
sum elastic_price_m1,d
bysort country year: asgen avgela_price_m1=elastic_price_m1, weight(reg)
sum avgela_price_m1,d	



* 1.2 model2:two-level nest: segment-->origin, ivfe, model-bodytype-trimlevel FE, segment-year
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize    ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m2=i.country2#i.year fe_mb_m2=i.groupmodel2#i.bodytype_harmonized4#i.trimlevel2    ///
	 fe_sy_m2=i.segment2#i.year, resid(resid_m2)) first ffirst cluster(group_mdltrim)
est store m2
matrix beta_m2=e(b)'
matrix list beta_m2

gen gamma_jcy_m2=logshare_jcy-beta_m2[1,1]*price-beta_m2[2,1]*logshare_cyjinso-beta_m2[3,1]*logshare_cyoins   ///
	-beta_m2[4,1]*tax-beta_m2[5,1]*energycost1-beta_m2[6,1]*loghw-beta_m2[7,1]*logw-fe_cy_m2
gen wtp_jcy_m2=gamma_jcy_m2/beta_m2[1,1]*-1

//elasticity and wtp for attributes
gen elastic_price_m2=((1-share_cyjinso)/(1-beta_m2[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m2[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m2[1,1]*price
sum elastic_price_m2,d
bysort country year: asgen avgela_price_m2=elastic_price_m2, weight(reg)
sum avgela_price_m2,d	


* 1.3 model 3: two-level nest: segment-->origin, ivfe, model-bodytype-trimlevel-segment FE, segment-year
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m3=i.country2#i.year fe_mb_m3=i.groupmodel2#i.bodytype_harmonized4#i.trimlevel2#i.segment2   ///
	 fe_sy_m3=i.segment2#i.year, resid(resid_m3)) first ffirst cluster(group_mdltrim)
est store m3
matrix beta_m3=e(b)'
matrix list beta_m3

gen gamma_jcy_m3=logshare_jcy-beta_m3[1,1]*price-beta_m3[2,1]*logshare_cyjinso-beta_m3[3,1]*logshare_cyoins    ///
	-beta_m3[4,1]*tax-beta_m3[5,1]*energycost1-beta_m3[6,1]*loghw-beta_m3[7,1]*logw-fe_cy_m3
gen wtp_jcy_m3=gamma_jcy_m3/beta_m3[1,1]*-1

//elasticity and wtp for attributes
gen elastic_price_m3=((1-share_cyjinso)/(1-beta_m3[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m3[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m3[1,1]*price
sum elastic_price_m3,d
bysort country year: asgen avgela_price_m3=elastic_price_m3, weight(reg)
sum avgela_price_m3,d	


* 1.4 model 4: two-level nest: segment-->origin, ivfe, model-bodytype-trimlevel-segment-fueltype FE, segment-year
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m4=i.country2#i.year fe_mb_m4=i.groupmodel2#i.bodytype_harmonized4#i.trimlevel2#i.segment2#i.fuelcat_harmonized2   ///
	 fe_sy_m4=i.segment2#i.year, resid(resid_m4)) first ffirst cluster(group_mdltrim)
est store m4
matrix beta_m4=e(b)'
matrix list beta_m4

gen gamma_jcy_m4=logshare_jcy-beta_m4[1,1]*price-beta_m4[2,1]*logshare_cyjinso-beta_m4[3,1]*logshare_cyoins    ///
	-beta_m4[4,1]*tax-beta_m4[5,1]*energycost1-beta_m4[6,1]*loghw-beta_m4[7,1]*logw-fe_cy_m4
gen wtp_jcy_m4=gamma_jcy_m4/beta_m4[1,1]*-1

//elasticity and wtp for attributes
gen elastic_price_m4=((1-share_cyjinso)/(1-beta_m4[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m4[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m4[1,1]*price
sum elastic_price_m4,d
bysort country year: asgen avgela_price_m4=elastic_price_m4, weight(reg)
sum avgela_price_m4,d	

* 1.5 model 5: two-level nest: segment-->origin, ivfe, vehicle FE, segment-year
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m5=i.country2#i.year fe_mb_m5=vehicle   ///
	 fe_sy_m5=i.segment2#i.year, resid(resid_m5)) first ffirst cluster(group_mdltrim)
est store m5
matrix beta_m5=e(b)'
matrix list beta_m5

gen gamma_jcy_m5=logshare_jcy-beta_m5[1,1]*price-beta_m5[2,1]*logshare_cyjinso-beta_m5[3,1]*logshare_cyoins    ///
	-beta_m5[4,1]*tax-beta_m5[5,1]*energycost1-beta_m5[6,1]*loghw-beta_m5[7,1]*logw-fe_cy_m5
gen wtp_jcy_m5=gamma_jcy_m5/beta_m5[1,1]*-1

//elasticity and wtp for attributes
gen elastic_price_m5=((1-share_cyjinso)/(1-beta_m5[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m5[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m5[1,1]*price
sum elastic_price_m5,d
bysort country year: asgen avgela_price_m5=elastic_price_m5, weight(reg)
sum avgela_price_m5,d	

* 1.6 model 6:model-bodytype FE, segment-year, groupmodel2-year, cy
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m6=i.country2#i.year fe_mb_m6=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_m6=i.segment2#i.year fe_my6=i.groupmodel2#i.year, resid(resid_m6)) first ffirst cluster(group_mdltrim)
est store m6
matrix beta_m6=e(b)'
matrix list beta_m6

gen gamma_jcy_m6=logshare_jcy-beta_m6[1,1]*price-beta_m6[2,1]*logshare_cyjinso-beta_m6[3,1]*logshare_cyoins    ///
	-beta_m6[4,1]*tax-beta_m6[5,1]*energycost1-beta_m6[6,1]*loghw-beta_m6[7,1]*logw-fe_cy_m6
gen wtp_jcy_m6=gamma_jcy_m6/beta_m6[1,1]*-1

//elasticity and wtp for attributes
gen elastic_price_m6=((1-share_cyjinso)/(1-beta_m6[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m6[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m6[1,1]*price
sum elastic_price_m6,d
bysort country year: asgen avgela_price_m6=elastic_price_m6, weight(reg)
sum avgela_price_m6,d	

* 1.7 model 7:model-bodytype FE, segment-year, bodytype-year, cy
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m7=i.country2#i.year fe_mb_m7=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_m7=i.segment2#i.year fe_by7=i.bodytype_harmonized4#i.year, resid(resid_m7)) first ffirst cluster(group_mdltrim)
est store m7
matrix beta_m7=e(b)'
matrix list beta_m7

gen gamma_jcy_m7=logshare_jcy-beta_m7[1,1]*price-beta_m7[2,1]*logshare_cyjinso-beta_m7[3,1]*logshare_cyoins    ///
	-beta_m7[4,1]*tax-beta_m7[5,1]*energycost1-beta_m7[6,1]*loghw-beta_m7[7,1]*logw-fe_cy_m7
gen wtp_jcy_m7=gamma_jcy_m7/beta_m7[1,1]*-1

//elasticity and wtp for attributes
gen elastic_price_m7=((1-share_cyjinso)/(1-beta_m7[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m7[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m7[1,1]*price
sum elastic_price_m7,d
bysort country year: asgen avgela_price_m7=elastic_price_m7, weight(reg)
sum avgela_price_m7,d	

* 1.8 model 8:model-bodytype FE, segment-year, trimlevel-year, cy
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m8=i.country2#i.year fe_mb_m8=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_m8=i.segment2#i.year fe_my8=i.trimlevel2#i.year, resid(resid_m8)) first ffirst cluster(group_mdltrim)
est store m8
matrix beta_m8=e(b)'
matrix list beta_m8

gen gamma_jcy_m8=logshare_jcy-beta_m8[1,1]*price-beta_m8[2,1]*logshare_cyjinso-beta_m8[3,1]*logshare_cyoins    ///
	-beta_m8[4,1]*tax-beta_m8[5,1]*energycost1-beta_m8[6,1]*loghw-beta_m8[7,1]*logw-fe_cy_m8
gen wtp_jcy_m8=gamma_jcy_m8/beta_m8[1,1]*-1

//elasticity and wtp for attributes
gen elastic_price_m8=((1-share_cyjinso)/(1-beta_m8[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m8[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m8[1,1]*price
sum elastic_price_m8,d
bysort country year: asgen avgela_price_m8=elastic_price_m8, weight(reg)
sum avgela_price_m8,d	

* 1.9 model 9:model-bodytype FE, segment-year, fuelcat-year, cy
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m9=i.country2#i.year fe_mb_m9=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_m9=i.segment2#i.year fe_fy9=i.fuelcat_harmonized2#i.year, resid(resid_m9)) first ffirst cluster(group_mdltrim)
est store m9
matrix beta_m9=e(b)'
matrix list beta_m9

gen gamma_jcy_m9=logshare_jcy-beta_m9[1,1]*price-beta_m9[2,1]*logshare_cyjinso-beta_m9[3,1]*logshare_cyoins    ///
	-beta_m9[4,1]*tax-beta_m9[5,1]*energycost1-beta_m9[6,1]*loghw-beta_m9[7,1]*logw-fe_cy_m9
gen wtp_jcy_m9=gamma_jcy_m9/beta_m9[1,1]*-1

//elasticity and wtp for attributes
gen elastic_price_m9=((1-share_cyjinso)/(1-beta_m9[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m9[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m9[1,1]*price
sum elastic_price_m9,d
bysort country year: asgen avgela_price_m9=elastic_price_m9, weight(reg)
sum avgela_price_m9,d	


* 1.10 model 10:model-bodytype FE, segment-year, groupmodel2-bodytype-year, cy
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m10=i.country2#i.year fe_mb_m10=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_m10=i.segment2#i.year fe_mby10=i.groupmodel2#i.bodytype_harmonized4#i.year, resid(resid_m10)) first ffirst cluster(group_mdltrim)
est store m10
matrix beta_m10=e(b)'
matrix list beta_m10

gen gamma_jcy_m10=logshare_jcy-beta_m10[1,1]*price-beta_m10[2,1]*logshare_cyjinso-beta_m10[3,1]*logshare_cyoins    ///
	-beta_m10[4,1]*tax-beta_m10[5,1]*energycost1-beta_m10[6,1]*loghw-beta_m10[7,1]*logw-fe_cy_m10
gen wtp_jcy_m10=gamma_jcy_m10/beta_m10[1,1]*-1

//elasticity and wtp for attributes
gen elastic_price_m10=((1-share_cyjinso)/(1-beta_m10[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m10[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m10[1,1]*price
sum elastic_price_m10,d
bysort country year: asgen avgela_price_m10=elastic_price_m10, weight(reg)
sum avgela_price_m10,d	

* 1.11 model 11:model-bodytype FE, segment-year, groupmodel2-year, bodytype-year, cy
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m11=i.country2#i.year fe_mb_m11=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_m11=i.segment2#i.year fe_my11=i.groupmodel2#i.year fe_by11=i.bodytype_harmonized4#i.year, resid(resid_m11)) first ffirst cluster(group_mdltrim)
est store m11
matrix beta_m11=e(b)'
matrix list beta_m11

gen gamma_jcy_m11=logshare_jcy-beta_m11[1,1]*price-beta_m11[2,1]*logshare_cyjinso-beta_m11[3,1]*logshare_cyoins    ///
	-beta_m11[4,1]*tax-beta_m11[5,1]*energycost1-beta_m11[6,1]*loghw-beta_m11[7,1]*logw-fe_cy_m11
gen wtp_jcy_m11=gamma_jcy_m11/beta_m11[1,1]*-1

//elasticity and wtp for attributes
gen elastic_price_m11=((1-share_cyjinso)/(1-beta_m11[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m11[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m11[1,1]*price
sum elastic_price_m11,d
bysort country year: asgen avgela_price_m11=elastic_price_m11, weight(reg)
sum avgela_price_m11,d	

* 1.12 model 12:model-bodytype FE, segment-year, groupmodel2-year, bodytype-year, trimlevel-year, cy
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m12=i.country2#i.year fe_mb_m12=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_m12=i.segment2#i.year fe_my12=i.groupmodel2#i.year fe_by12=i.bodytype_harmonized4#i.year fe_ty12=i.trimlevel2#i.year, resid(resid_m12)) first ffirst cluster(group_mdltrim)
est store m12
matrix beta_m12=e(b)'
matrix list beta_m12

gen gamma_jcy_m12=logshare_jcy-beta_m12[1,1]*price-beta_m12[2,1]*logshare_cyjinso-beta_m12[3,1]*logshare_cyoins    ///
	-beta_m12[4,1]*tax-beta_m12[5,1]*energycost1-beta_m12[6,1]*loghw-beta_m12[7,1]*logw-fe_cy_m12
gen wtp_jcy_m12=gamma_jcy_m12/beta_m12[1,1]*-1

//elasticity and wtp for attributes
gen elastic_price_m12=((1-share_cyjinso)/(1-beta_m12[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m12[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m12[1,1]*price
sum elastic_price_m12,d
bysort country year: asgen avgela_price_m12=elastic_price_m12, weight(reg)
sum avgela_price_m12,d	

* 1.13 model 13:model-bodytype FE, segment-year, groupmodel2-year, bodytype-year, trimlevel-year, fuelcat-year, cy
sort country year vehicle
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(fe_cy_m13=i.country2#i.year fe_mb_m13=i.groupmodel2#i.bodytype_harmonized4    ///
	 fe_sy_m13=i.segment2#i.year fe_my13=i.groupmodel2#i.year fe_by13=i.bodytype_harmonized4#i.year fe_ty13=i.trimlevel2#i.year fe_fy13=i.fuelcat_harmonized2#i.year, resid(resid_m13)) first ffirst cluster(group_mdltrim)
est store m13
matrix beta_m13=e(b)'
matrix list beta_m13

gen gamma_jcy_m13=logshare_jcy-beta_m13[1,1]*price-beta_m13[2,1]*logshare_cyjinso-beta_m13[3,1]*logshare_cyoins    ///
	-beta_m13[4,1]*tax-beta_m13[5,1]*energycost1-beta_m13[6,1]*loghw-beta_m13[7,1]*logw-fe_cy_m13
gen wtp_jcy_m13=gamma_jcy_m13/beta_m13[1,1]*-1

//elasticity and wtp for attributes
gen elastic_price_m13=((1-share_cyjinso)/(1-beta_m13[2,1])+(1-share_cyoins)*share_cyjinso/(1-beta_m13[3,1])    ///
	+(1-share_cys)*share_cyoins*share_cyjinso)*beta_m13[1,1]*price
sum elastic_price_m13,d
bysort country year: asgen avgela_price_m13=elastic_price_m13, weight(reg)
sum avgela_price_m13,d	

save "demand_appendix_different_FEs_20210418.dta", replace


* 2. DID: using quality from different demand models with different FEs
use "demand_appendix_different_FEs_20210418.dta", clear
foreach x in m1 m2 m3 m4 m5{
	reghdfe wtp_jcy_`x' c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3   ///
		[fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year) vce(cluster group_mdltrim) 
	est store wtp_`x'
	count if e(sample)==1
	matrix n_wtp_`x'=r(N)
	test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
	matrix test_wtp_`x'=(r(drop),  r(df_r), r(F), r(df), r(p))'
}


//model 6: mb, sy, my, cy
reghdfe wtp_jcy_m6 c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3 [fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year i.groupmodel2#i.year) vce(cluster group_mdltrim) 
est store wtp_m6
count if e(sample)==1
matrix n_wtp_m6=r(N)
test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
matrix test_wtp_m6=(r(drop),  r(df_r), r(F), r(df), r(p))'

//model 7: mb, sy, bodytype-year, cy
reghdfe wtp_jcy_m7 c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3 [fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year i.bodytype_harmonized4#i.year) vce(cluster group_mdltrim) 
est store wtp_m7
count if e(sample)==1
matrix n_wtp_m7=r(N)
test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
matrix test_wtp_m7=(r(drop),  r(df_r), r(F), r(df), r(p))'

//model 8: mb, sy, trim-year, cy
reghdfe wtp_jcy_m8 c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3 [fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year i.trimlevel2#i.year) vce(cluster group_mdltrim) 
est store wtp_m8
count if e(sample)==1
matrix n_wtp_m8=r(N)
test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
matrix test_wtp_m8=(r(drop),  r(df_r), r(F), r(df), r(p))'

//model 9: mb, sy, fuelcat-year, cy
reghdfe wtp_jcy_m9 c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3 [fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year i.fuelcat_harmonized2#i.year) vce(cluster group_mdltrim) 
est store wtp_m9
count if e(sample)==1
matrix n_wtp_m9=r(N)
test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
matrix test_wtp_m9=(r(drop),  r(df_r), r(F), r(df), r(p))'

//model 10: mb, sy, model-body type-year, cy
reghdfe wtp_jcy_m10 c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3 [fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year i.groupmodel2#i.bodytype_harmonized4#i.year) vce(cluster group_mdltrim) 
est store wtp_m10
count if e(sample)==1
matrix n_wtp_m10=r(N)
test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
matrix test_wtp_m10=(r(drop),  r(df_r), r(F), r(df), r(p))'

//model 11: mb, sy, cy, my, by
reghdfe wtp_jcy_m11 c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3 [fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year i.groupmodel2#i.year i.bodytype_harmonized4#i.year) vce(cluster group_mdltrim) 
est store wtp_m11
count if e(sample)==1
matrix n_wtp_m11=r(N)
test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
matrix test_wtp_m11=(r(drop),  r(df_r), r(F), r(df), r(p))'

//model 12: mb, sy, cy, my, by, ty
reghdfe wtp_jcy_m12 c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3 [fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year i.groupmodel2#i.year i.bodytype_harmonized4#i.year i.trimlevel2#i.year) vce(cluster group_mdltrim) 
est store wtp_m12
count if e(sample)==1
matrix n_wtp_m12=r(N)
test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
matrix test_wtp_m12=(r(drop),  r(df_r), r(F), r(df), r(p))'

//model 13: mb, sy, cy, my, by, ty, fy
reghdfe wtp_jcy_m13 c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3 [fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year i.groupmodel2#i.year i.bodytype_harmonized4#i.year i.trimlevel2#i.year i.fuelcat_harmonized2#i.year) vce(cluster group_mdltrim) 
est store wtp_m13
count if e(sample)==1
matrix n_wtp_m13=r(N)
test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
matrix test_wtp_m13=(r(drop),  r(df_r), r(F), r(df), r(p))'

//output
esttab wtp_m1 wtp_m2 wtp_m3 wtp_m4 wtp_m5, b(%7.3f) se(%7.3f) stats(N r2_a) nogap br compress order(*.post#* stringency_* _cons) star(* 0.1 ** 0.05 *** 0.01)

esttab wtp_m6 wtp_m7 wtp_m8 wtp_m9 wtp_m10 wtp_m11 wtp_m12 wtp_m13, b(%7.3f) se(%7.3f) stats(N r2_a) nogap br compress order(*.post#* stringency_* _cons) star(* 0.1 ** 0.05 *** 0.01)

matrix drop A 
matrix A=(test_wtp_m1, test_wtp_m2, test_wtp_m3, test_wtp_m4, test_wtp_m5, test_wtp_m6, test_wtp_m7, test_wtp_m8, test_wtp_m9, test_wtp_m10, test_wtp_m11, test_wtp_m12, test_wtp_m13  )
matrix list A
matrix B=(n_wtp_m1, n_wtp_m2, n_wtp_m3, n_wtp_m4, n_wtp_m5, n_wtp_m6, n_wtp_m7, n_wtp_m8, n_wtp_m9, n_wtp_m10, n_wtp_m11, n_wtp_m12, n_wtp_m13 )
matrix list B

* 2.2 sy, vehicle, cy fixed effects
foreach x in m1 m2 m3 m4 m5 m6 m7 m8 m9 m10 m11 m12 m13{
	reghdfe wtp_jcy_`x' c.stringency_f_m3#2.post c.stringency_f_m3#3.post stringency_f_m3   ///
		[fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year) vce(cluster group_mdltrim) 
	est store wtp_`x'
	count if e(sample)==1
	matrix n_wtp_`x'=r(N)
	test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
	matrix test_wtp_`x'=(r(drop),  r(df_r), r(F), r(df), r(p))'
}

//output
esttab wtp_m1 wtp_m2 wtp_m3 wtp_m4 wtp_m5, b(%7.3f) se(%7.3f) stats(N r2_a) nogap br compress order(*.post#* stringency_* _cons) star(* 0.1 ** 0.05 *** 0.01)

esttab wtp_m6 wtp_m7 wtp_m8 wtp_m9 wtp_m10 wtp_m11 wtp_m12 wtp_m13, b(%7.3f) se(%7.3f) stats(N r2_a) nogap br compress order(*.post#* stringency_* _cons) star(* 0.1 ** 0.05 *** 0.01)

matrix drop A B
matrix A=(test_wtp_m1, test_wtp_m2, test_wtp_m3, test_wtp_m4, test_wtp_m5, test_wtp_m6, test_wtp_m7, test_wtp_m8, test_wtp_m9, test_wtp_m10, test_wtp_m11, test_wtp_m12, test_wtp_m13  )
matrix list A
matrix B=(n_wtp_m1, n_wtp_m2, n_wtp_m3, n_wtp_m4, n_wtp_m5, n_wtp_m6, n_wtp_m7, n_wtp_m8, n_wtp_m9, n_wtp_m10, n_wtp_m11, n_wtp_m12, n_wtp_m13 )
matrix list B



********************************************************************************
**# Appendix Table 28: Quality regression, Re-design
********************************************************************************
use "data.dta", clear
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(i.country2#i.year i.groupmodel2#i.bodytype_harmonized4    ///
	 i.segment2#i.year) first ffirst cluster(group_mdltrim)
est store nl
matrix beta=e(b)'
matrix list beta
gen snl=e(sample)

//normalize the loghw and logw
drop loghw3 logw2
gen loghw3=loghw/abs(beta[1,1])
gen logw2=logw/abs(beta[1,1])
label var loghw3 "loghw (hp/kg) normalized by the absolute value of price coefficient"
label var logw2 "logw (ton) normalized by the absolute value of price coefficient"

isid vehicle country year
* 1. compute the percentage change in each year--vehicle level
sum fuelcons enginehorsepower size grossvehicleweight if snl==1
sum fuelcons enginehorsepower size grossvehicleweight

foreach x in fuelcons enginehorsepower size grossvehicleweight{
	bysort vehicle year: egen avg_`x'=mean(`x')			//unweighted average
}

duplicates drop vehicle year, force

bysort vehicle: egen appear_year=min(year)
bysort vehicle: egen exit_year=max(year)
bysort vehicle: gen duration=exit_year-appear_year+1
tab appear_year
tab exit_year
tab duration

sort vehicle year
foreach x in fuelcons enginehorsepower size grossvehicleweight{
	sort vehicle year
	bysort vehicle: gen change0_`x'=avg_`x'-avg_`x'[_n-1]
	sort vehicle year
	bysort vehicle: gen pc0_`x'=(avg_`x'-avg_`x'[_n-1])/avg_`x'[_n-1]
	gen change_`x'=abs(change0_`x')
	gen pc_`x'=abs(pc0_`x')
	replace change0_`x'=0 if year==appear_year 
	replace pc0_`x'=0 if year==appear_year 
	replace change_`x'=0 if year==appear_year 
	replace pc_`x'=0 if year==appear_year 	
}	

sum change_* pc_* change0_* pc0_*
sort vehicle year 
br vehicle year avg_fuelcons change0_fuelcons change_fuelcons pc0_fuelcons pc_fuelcons
br vehicle year avg_fuelcons change0_fuelcons change_fuelcons pc0_fuelcons pc_fuelcons if pc_fuelcons>=.

sort vehicle year
foreach x in fuelcons enginehorsepower size grossvehicleweight{
	sum pc_`x', d
	gen flgpc_`x'=1 if pc_`x'>r(p99) & pc_`x'<.
	replace flgpc_`x'=0 if flgpc_`x'==.
}
sum pc_fuelcons if flgpc_fuelcons==0, d
matrix A=(r(N), r(mean), r(sd), r(min), r(max), r(p1),  r(p5),  r(p10), r(p25), r(p50), r(p75), r(p90), r(p95), r(p99))
matrix list A
foreach x in enginehorsepower size grossvehicleweight{
	sum pc_`x' if flgpc_`x'==0,d
	matrix A=(A\(r(N), r(mean), r(sd), r(min), r(max), r(p1),  r(p5),  r(p10), r(p25), r(p50), r(p75), r(p90), r(p95), r(p99)))
}
matrix list A

sort vehicle year
br vehicle year avg_fuelcons pc_fuelcons if (pc_fuelcons>A[1,14]) & pc_fuelcons<.
twoway hist pc_fuelcons if flgpc_fuelcons==0, fcolor(none) 


* 2. generate the "re-design" variable: vehicle-level 
* Note: if a vehicle is re-designed, then one of the four variables changes more than x percentage
//fuelcons 
gen redesign90_fuelcons=cond(pc_fuelcons>A[1,12] & pc_fuelcons<. & flgpc_fuelcons==0,1,0)
//hp 
gen redesign90_enginehorsepower=cond(pc_enginehorsepower>A[2,12] & pc_enginehorsepower<. & flgpc_enginehorsepower==0,1,0)
//size
gen redesign90_size=cond(pc_size>A[3,12] & pc_size<. & flgpc_size==0,1,0)
//weight
gen redesign90_grossvehicleweight=cond(pc_grossvehicleweight>A[4,12] & pc_grossvehicleweight<. & flgpc_grossvehicleweight==0,1,0)
//"re-design" variables
egen redesignnum90=anycount(redesign90_fuelcons redesign90_enginehorsepower redesign90_size redesign90_grossvehicleweight), values(1)
egen redesign90=rowmax(redesign90_fuelcons redesign90_enginehorsepower redesign90_size redesign90_grossvehicleweight)
sum redesign*

save "vehicleyear_redesign" , replace
duplicates drop groupmodel2 year, force 

* 3. compute the percentage change in each year--model level, method 1
//take the average of the percentage change within a model 
foreach x in fuelcons enginehorsepower size grossvehicleweight{
	bysort groupmodel2 year: egen changemm_`x'=mean(change_`x')
	bysort groupmodel2 year: egen pcmm_`x'=mean(pc_`x')
}
sort groupmodel2 year 
sum pcmm_fuelcons, d
matrix A=(r(N), r(mean), r(sd), r(min), r(max), r(p1),  r(p5),  r(p10), r(p25), r(p50), r(p75), r(p90), r(p95), r(p99))
matrix list A
foreach x in enginehorsepower size grossvehicleweight{
	sum pcmm_`x' ,d
	matrix A=(A\(r(N), r(mean), r(sd), r(min), r(max), r(p1),  r(p5),  r(p10), r(p25), r(p50), r(p75), r(p90), r(p95), r(p99)))
}
matrix list A

* 4. generate the "re-design" variable: model-level
//fuelcons 
gen redesignmm90_fuelcons=cond(pcmm_fuelcons>A[1,12] & pcmm_fuelcons<.,1,0)
//hp 
gen redesignmm90_enginehorsepower=cond(pcmm_enginehorsepower>A[2,12] & pcmm_enginehorsepower<.,1,0)
//size
gen redesignmm90_size=cond(pcmm_size>A[3,12] & pcmm_size<.,1,0)
//weight 
gen redesignmm90_grossvehicleweight=cond(pcmm_grossvehicleweight>A[4,12] & pcmm_grossvehicleweight<.,1,0)
//"re-design" variables
egen redesignmmnum90=anycount(redesignmm90_fuelcons redesignmm90_enginehorsepower redesignmm90_size redesignmm90_grossvehicleweight), values(1)

egen redesignmm90=rowmax(redesignmm90_fuelcons redesignmm90_enginehorsepower redesignmm90_size redesignmm90_grossvehicleweight)

sum redesignmm*

* 5. drop extreme values 
sort groupmodel2 year
foreach x in fuelcons enginehorsepower size grossvehicleweight{
	sum pcmm_`x', d
	gen flgpcmm_`x'=1 if pcmm_`x'>r(p99) & pcmm_`x'<.
	replace flgpcmm_`x'=0 if flgpcmm_`x'==.
}
sum pcmm_fuelcons if flgpcmm_fuelcons==0, d
matrix A=(r(N), r(mean), r(sd), r(min), r(max), r(p1),  r(p5),  r(p10), r(p25), r(p50), r(p75), r(p90), r(p95), r(p99))
matrix list A
foreach x in enginehorsepower size grossvehicleweight{
	sum pcmm_`x' if flgpcmm_`x'==0,d
	matrix A=(A\(r(N), r(mean), r(sd), r(min), r(max), r(p1),  r(p5),  r(p10), r(p25), r(p50), r(p75), r(p90), r(p95), r(p99)))
}
matrix list A

//fuelcons 
gen redesignmm290_fuelcons=cond(pcmm_fuelcons>A[1,12] & pcmm_fuelcons<. & flgpcmm_fuelcons==0,1,0)
//hp 
gen redesignmm290_enginehorsepower=cond(pcmm_enginehorsepower>A[2,12] & pcmm_enginehorsepower<. & flgpcmm_enginehorsepower==0,1,0)
//size
gen redesignmm290_size=cond(pcmm_size>A[3,12] & pcmm_size<. & flgpcmm_size==0,1,0)
//weight 
gen redesignmm290_grossvehicleweight=cond(pcmm_grossvehicleweight>A[4,12] & pcmm_grossvehicleweight<.  & flgpcmm_grossvehicleweight==0,1,0)
//"re-design" variables
egen redesignmm2num90=anycount(redesignmm290_fuelcons redesignmm290_enginehorsepower redesignmm290_size redesignmm290_grossvehicleweight), values(1)

egen redesignmm290=rowmax(redesignmm290_fuelcons redesignmm290_enginehorsepower redesignmm290_size redesignmm290_grossvehicleweight)

sum redesignmm*
save "groupmodelyear_redesign1" , replace


* 6. DDD with re-design variables : dropping extreme values
//merge re-design variables with the original data
use "data.dta", clear
ivreghdfe logshare_jcy tax energycost1 loghw logw logsize  ///
	(price logshare_cyjinso logshare_cyoins=sum_length_same sum_length_diff sum_width_same sum_width_diff sum_height_same sum_height_diff    ///
	sum_enginecyl_same sum_enginecyl_diff    ///
	sum_length_oins1 sum_width_oins1 sum_height_oins1    ///
	sum_length_oins2 sum_width_oins2 sum_height_oins2    ///
	sum_enginecyl_oins1 sum_enginecyl_oins2),   ///
	absorb(i.country2#i.year i.groupmodel2#i.bodytype_harmonized4    ///
	 i.segment2#i.year) first ffirst cluster(group_mdltrim)
est store nl
matrix beta=e(b)'
matrix list beta
gen snl=e(sample)

//normalize the loghw and logw
drop loghw3 logw2
gen loghw3=loghw/abs(beta[1,1])
gen logw2=logw/abs(beta[1,1])
label var loghw3 "loghw (hp/kg) normalized by the absolute value of price coefficient"
label var logw2 "logw (ton) normalized by the absolute value of price coefficient"

//merge re-design informations 
merge m:1 vehicle year using "vehicleyear_redesign", keepusing(redesign?? redesignnum?? flgpc_*)
drop _merge
sum redesign?? redesignnum?? flgpc_*
merge m:1 groupmodel2 year using "groupmodelyear_redesign1", keepusing(redesignmm* flgpcmm*)
drop _merge
sum redesignmm* flgpcmm*
sort vehicle country year

// DDD: model-level, method 1, drop extreme values
foreach y in 90{
	foreach x in wtp_jcy1 loghw3 logw2 logprice {
		reghdfe `x' c.stringency_f_m3##i.post##i.redesignmm2`y' [fw=reg], absorb(vehicle i.country2 i.year i.country2#i.year i.segment2#i.year) vce(cluster group_mdltrim) 
		est store wf0`y'_`x'
		count if e(sample)==1
		matrix n_wf0`y'_`x'=r(N)
		test c.stringency_f_m3#2.post c.stringency_f_m3#3.post
		matrix test_wf0`y'_`x'=(r(drop),  r(df_r), r(F), r(df), r(p))'
		test c.stringency_f_m3#2.post#1.redesignmm2`y' c.stringency_f_m3#3.post#1.redesignmm2`y'
		matrix test2_wf0`y'_`x'=(r(drop),  r(df_r), r(F), r(df), r(p))'
	}
}

//90%
esttab wf090_wtp_jcy1 wf090_loghw3 wf090_logw2 wf090_logprice,     ///
	b(%7.3f) se(%7.3f) stats(N r2_a) br compress order(2.post#1.redesignmm2??#c.stringency_f_m3 3.post#1.redesignmm2??#c.stringency_f_m3 2.post#c.stringency_f_m3  3.post#c.stringency_f_m3 2.post#1.redesignmm2?? 3.post#1.redesignmm2?? 1.redesignmm2??#c.stringency_f_m3 1.redesignmm2??  stringency_* _cons) keep(2.post#1.redesignmm2??#c.stringency_f_m3 3.post#1.redesignmm2??#c.stringency_f_m3 2.post#c.stringency_f_m3  3.post#c.stringency_f_m3 2.post#1.redesignmm2?? 3.post#1.redesignmm2?? 1.redesignmm2??#c.stringency_f_m3 1.redesignmm2?? ) star(* 0.1 ** 0.05 *** 0.01) nogap noline

matrix A=(test_wf090_wtp_jcy1, test_wf090_loghw3, test_wf090_logw2, test_wf090_logprice)
matrix C=(test2_wf090_wtp_jcy1, test2_wf090_loghw3, test2_wf090_logw2, test2_wf090_logprice)
matrix B=(n_wf090_wtp_jcy1, n_wf090_loghw3, n_wf090_logw2, n_wf090_logprice)
matrix D=(B\C[3,1...]\ C[5,1...]\A[3,1...]\A[5,1...])
matrix list D
matrix drop A B C D

































